Classical Algorithm for Generating Multi-Mode Bosonic Transition Spectra by Phase Space Sampling

ABSTRACT

This disclosure relates to a method and system for efficiently and analytically generating spectra associated with transitions between thermal states in a multi-mode bosonic system using a classical algorithm to compute Fourier components of the transition spectra of the multi-mode bosonic system based on a representation space sampling method inspired by quantum Gaussian boson sampling followed by an inverse Fourier transform. The disclosed method and system particularly apply to efficient estimation of molecular vibronic spectra. For example, an exact solution of the Fourier components of molecular vibronic spectra at zero temperature may be analytically computed in a representation space by using, for example, a positive P-representation of the multi-mode bosonic vibronic quantum states of the molecule. Such a method may be further applied to more general vibronic spectroscopy, such as computing molecular vibronic spectra at finite temperatures by introducing additional auxiliary bosonic modes and computing vibronic spectra associated with non-thermal states, such as Fock states.

CROSS REFERENCE

This application is based on and claims the benefit of priority to U.S.Provisional Patent Application Nos. 63/300,320 and 63/393,048 filed onJan. 18, 2022 and Jul. 28, 2022, respectively. These prior provisionalpatent applications are herein incorporated by reference in theirentireties.

GOVERNMENT FUNDING

This invention was made with government support under grant numbers W911NF-18-1-0020, W911 NF-18-1-0212, W911 NF-16-1-0349, and W911NF-21-1-0325 awarded by the Army Research Office, grant numbersFA9550-19-1-0399, FA9550-21-1-0209, FA8649-21-P-0781, FA9550-18-1-0148,and FA9550-21-1-0008 awarded by Air Force Office of Scientific Research,grant numbers 1640959, 1936118, 1941583, 2137642, and 2044923 awarded bythe National Science Foundation, and grant number DE-SC0020360 awardedby the Department of Energy. The government has certain rights in theinvention.

FIELD OF THE INVENTION

This disclosure relates to a method and system for efficientlygenerating spectra associated with transitions between thermal states ina multi-mode bosonic system using a classical algorithm based on a phasespace sampling method inspired by a quantum Gaussian boson samplingmethodology. The disclosed method and system particularly apply toefficient estimation of molecular vibronic spectra.

BACKGROUND OF THE INVENTION

Vibronic spectra represent an essential and characteristic property of amolecule composition and structure. As such, molecular vibronic spectramay be relied upon to analyze chemical composition of molecules and tostudy molecular structures. Molecular vibronic spectra of a plurality ofreference molecular structures may be computationally generated and usedas a library of molecular fingerprints for experimentally identifyingand studying molecular compounds. However, the computation of vibronicspectra of a particular molecular structure has been known to be achallenging task, with its complexity rapidly growing with the sizeand/or the number of vibrational modes of the molecular structure. Forexample, the best-known classical algorithm scales combinatorially withthe size of the molecular system. While a quantum simulator,particularly the ones based on Gaussian boson sampling of multi-modephotons, may be implemented to efficiently generate the characteristicmolecular vibronic spectra, such a simulator does require a quantumcomputing machine, such as a linear optical circuits constructed toprocess delicate non-classical photonic states generated bynon-classical photon sources.

SUMMARY OF THE INVENTION

This disclosure relates to a method and system for efficiently andanalytically generating spectra associated with transitions betweenthermal states in a multi-mode bosonic system using a classicalalgorithm to compute Fourier components of the transition spectra of themulti-mode bosonic system based on a representation space samplingmethod inspired by quantum Gaussian boson sampling followed by aninverse Fourier transform. The disclosed method and system particularlyapply to efficient estimation of molecular vibronic spectra. Forexample, an exact solution of the Fourier components of molecularvibronic spectra at zero temperature may be analytically computed in aphase space by using, for example, a positive P-representation of themulti-mode bosonic vibronic quantum states of the molecule. Such amethod may be further applied to more general vibronic spectroscopy,such as computing molecular vibronic spectra at finite temperatures byintroducing additional auxiliary bosonic modes and computing vibronicspectra associated with non-thermal states.

In some example implementations, a method for estimating a transitionspectra in a bosonic system is disclosed. The bosonic system may beassociated with M bosonic modes, M being an integer equal to or largerthan 2. The method may include determining a unitary operation modifiedfrom a transition operator associated with the bosonic system;processing each of N sampling bosonic states of the bosonic system in arepresentation space of the bosonic system to generate N sets of samplesusing at least the unitary operation, wherein each set of samplescorrespond to one of the N sampling bosonic states and are sampled fromat least one variable in the representation space, N being a positiveinteger; analytically generating Fourier components the transitionspectra based on the N sets of samples; and inverse-transforming theFourier components to estimate the transition spectra.

In any one of the example implementations above, the unitary operationmay include a complex mixed position and momentum operator.

In any one of the example implementations above, each of the N samplingbosonic states may be a non-Gaussian state of the M bosonic modes.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the invention, reference is made tothe following description and accompanying drawings, in which:

FIG. 1 illustrates molecular vibronic spectra generation process basedon a quantum Gaussian algorithm and a quantum-inspired classicalalgorithm;

FIG. 2 shows an example data and logic flow for the quantum-inspiredclassical algorithm of FIG. 1 for generating molecular vibronic spectra;

FIG. 3 shows comparison of molecular vibronic spectra for CH₂O₂, 1¹A′→1²A′ as generated based on direct computation (ideal bars) and based onthe quantum-inspired classical and analytical algorithm of FIG. 1(“Analytical Solution” points and curve);

FIG. 4 illustrates transformation of molecular vibronic spectra atfinite temperature (upper) and treating the problem by introducingauxiliary modes and two-mode squeezing operation and using thezero-temperature algorithm (lower);

FIG. 5 illustrates transformation of molecular vibronic spectra with aFock state as an input;

FIG. 6 illustrates comparison of molecular vibronic spectra for CH₂O₂,1¹ A′→1² A′ as generated based on direct computation (ideal bars) andbased on Gaussian boson sampling simulation (“GBS” points and curve);

FIG. 7 illustrates comparison of total variance distance (TVD) betweenideal direct computation (upper curve with circular points) and thequantum-inspired classical algorithm (lower curve with triangularpoints) of FIG. 3 (the insert is in logarithm scale);

FIG. 8 illustrates virtual ideal vibronic spectra and empirical vibronicspectra generated using the classical algorithm with N=10⁵ samples andwith a total variation distance is of 0.00222;

FIG. 9 illustrates virtual ideal vibronic spectra and empirical vibronicspectra generated by Gaussian boson sampling experimental raw data withN=10⁶ samples and with a total variation distance of 0.01215; and

FIG. 10 illustrates Virtual FCP with a randomly generated weights forall 144 vibronic modes of a molecular system. The unit of wave number iscat. The generation of each instance using classical algorithm costsabout 90 secs using a single core computing system.

DETAILED DESCRIPTION Bosonic Systems and Molecular Vibronic Spectra

Atomic/nuclear components of a compound molecular structure interactwith one another and thus may vibrate in multiple number (e.g., integerM) of vibrational modes. Each vibrational mode, approximately, may bemodeled as a harmonic oscillator of a particular frequency. The multiplevibrational modes thus may correspond to a multitude of characteristicvibrational frequencies cv, within the molecular structure, where i∈[1,2, . . . , M]. Quantum mechanically, vibrations in each of the modes maybe described by quantum states |m_(i)

(where m_(i)∈[0, 1,]) with ladder-like energy levels each associatedwith a quantum number m_(i) and an energy determined by

ω_(i)(m_(i)−½) (for m_(i)=0, 1, 2, . . . , ∞), where the vacuum level isrepresented by

ω_(i)/2. An arbitrary molecular vibrational state thus may be describedby a quantum superposition of these multi-mode base states. Themolecular vibrations may thus be considered a bosonic system, similar toa system of multi-mode photons. For convenience, a molecular vibrationalquanta may be referred to as a photon in this disclosure.

The various vibrational modes are essential characteristics of amolecular structure. Molecular vibrations may be further coupled to theelectronic states or transitions in the molecule. For example,electronic transitions may change or perturb the nuclear structure ofthe molecule, thereby introducing new sets of vibrational modes. Thevibrational modes of a molecular may be investigated spectroscopically,for example, as fine spectral lines or structures within an electronictransition spectra. Such spectra may be referred to as molecularvibronic spectra. Vibronic spectra thus represent a fundamental propertyof a molecule. As such, molecular vibronic spectra may be relied upon toanalyze chemical compositions and study molecular structures. In manycircumstances, computation of such spectra may be critical in order toderive a molecular structure based on measured vibronic spectra. Forexample, molecular vibronic spectra of a plurality of referencemolecular structures may be computationally generated and used as alibrary of molecular fingerprints for experimentally identifying andstudying molecular compounds.

However, traditional classical computation of vibronic spectra of aparticular molecular structure has been known to be a challenging task,with its complexity rapidly growing with the size and/or the number ofthe vibrational modes of the molecular structure. For example, as shownbelow, the best-known classical algorithm scales combinatorially withthe size of the molecular system. While a quantum simulator,particularly the ones based on Gaussian boson sampling of multi-modephotons, may be implemented to efficiently generate the characteristicmolecular vibronic spectra, such a simulator does require a quantumcomputing machine, such as a linear optical circuits constructed toprocess delicate non-classical photonic states generated from anon-classical multi-mode photon source, which is usually challengingfrom experimental standpoint.

In this disclosure, a method and system for efficiently generatingmolecular vibrational spectra associated with transitions betweenvibrational thermal states in a molecular structure are described. Thedisclosed method and system are based on a quantum inspired classicalalgorithm. In particular, the disclosed classical algorithm analyticallyobtain Fourier components of the molecular vibrational spectra at zerotemperature by sampling phase-space variables in a positiveP-representation of a thermal multi-mode molecular vibrational state.Such a sampling technique is inspired by a quantum approach involvingGaussian boson sampling. While this algorithm relies on aquantum-inspired simulation, it nevertheless provides an exact andefficient solution of the Fourier components of molecular vibronicspectra in a classical manner without needing a quantum simulator. Themolecular vibrational spectra may then be generated by performinginverse Fourier transform of the Fourier components. Such a method isfurther extended to generate molecular vibronic spectra at finitetemperatures by introducing additional auxiliary modes and two-modesqueezing operations.

In more general cases, a Fock state rather than a thermal state as aninitial state may be considered. The Fock state may correspond to singlevibronic levels and may be written as a loop halfnian of a matrix. AMonte Carlo method using the positive P-representation may be used toestimate the molecular vibronic spectra in that case. However, the MonteCarlo method becomes inefficient when the number of quanta in theinitial Fock states grows. In some other implementations of a similarmethod setup in single-photon boson rather than multi-photon sampling toobtain corresponding virtual molecular vibronic spectra, thecorresponding Fourier components may become hard to compute (#P-hard),in stark contrast to Gaussian boson sampling which has an exactsolution, rendering a computational advantage for Gaussian bosonsampling over single-photon boson sampling. Nonetheless, a practicallyrelevant question is whether the Fourier components can be additivelyapproximated to within the same 1/poly(N) accuracy that can be obtainedusing N number of samples for a boson sampling experiment.Interestingly, when squeezing parameters are zero, the quantity reducesto a permanent, which can be approximated by Gurvits's algorithm withinan additive error. As such, an approximation method may be implementedto tackle the problem using the positive P-representation withpost-selection and another method that approximates a loop hafnian.While there may be a limitation of the method of positiveP-representation with post-selection as the number of quanta in theinitial Fock state grows, a regime may be identified where theapproximation error can be efficiently controlled even for large systemsize.

While the example implementations below are provided in the context ofmolecular vibronic spectral generation, the underlying principles andprocedure may be applicable to generating transition spectra in othermulti-mode bosonic systems.

Example Implementations of Quantum-Inspired Classical Algorithm forEstimating Molecular Vibronic Spectra

An example procedure for analytically generating molecular vibronicspectra using an efficient classical algorithm is shown in FIG. 1 incomparison to a typical quantum simulator for the same purposes. Inparticular, the analytical and classical algorithm is shown by the lowerdata-flow path 110 whereas the quantum simulation procedure isillustrated by the upper data-flow path 120. For both approaches, thecomputational output includes the vibronic spectra 130 for the zerotemperature molecule (e.g., absorption vibronic spectra), although theexample method/system disclosed herein can be expanded in non-zerotemperatures, as described in further detail below.

For the quantum simulation data path 120, an example Gaussian bosonsampling procedure may be performed. For example, a product quantumstate of single-mode Gaussian quantum states of multi-mode (e.g.,M-mode) photons may be initially determined, prepared, and generatedfrom appropriately configured and controlled photon sources. Such aproduct quantum state may represents a dressed and displaced initialstate (e.g., vacuum state for zero temperature and generated as an inputdressed and displaced M-mode photon state, as shown by 122, to a quantumcircuit 124. The quantum circuit 124 may include a plurality of beamsplitters and other linear optical components to simulate a unitaryquantum operation. A quantum unitary operator associated with thequantum circuit 124 may be determined in accordance with exampleimplementations described in further detail below.

For example, the unitary operator corresponding to the unitary quantumoperation of the quantum circuit 124, in combination with the dressingoperator and the displacement operator for generating the dressed anddisplaced initial quantum state of the input product quantum state ofthe M-mode photon state, may correspond to a transition operator of themolecular system between the ground state and vibronic thermal states.The linear optical components of the quantum processing circuit may thusbe configured to preform the operation corresponding to the unitaryoperator. The output of the linear operation may then be projected tothe various photon modes and detected as photon numbers in the M photonmodes, {m}. The array of linear optical components and the detection inthe M-photon modes are thus equivalent to a Gaussian Boson sampler, asshown by 124 and 126.

A set of N samples, represented by S={m^((i))}_(i=1) ^(N) and as shownby 128 in FIG. 1 , may be selected. The samples may then be postprocessed by first being sorted according to total energy of eachsample, where the total energy of each sample corresponds to the energyof all vibronic quanta (e.g., m^((i))·ω′, where ω′ represents the M-modefrequencies). In some practical implementations, for a particularspectra frequency, ω_(vib), a window Δω_(vib), may be introduced for thepost processing procedure. In other words, when a samples totalenergy-equivalent frequency m^((i))·ω′ lies between [ω_(vib)−Δω_(vib)/2,ω_(vib)+Δω_(vib)/2] for a certain ω_(vib), the count for a bincorresponding to ω_(vib), is incremented. The distribution of finalcounts in all ω_(vib) bins forms the simulated molecular vibronicspectra FCP(ω_(vib)), as long as sufficient N number of samples arecollected from the quantum linear optical circuit and are postprocessed.

In some example implementations, a classical algorithm may beimplemented as shown by the procedure 110 in FIG. 1 . While such analgorithm as described in further detail below is inspired by thequantum Gaussian boson sampling procedure described above (procedure 120of FIG. 1 ), it nevertheless would not need to rely on an actuallyquantum simulation machine. In order to generate molecular vibronicspectra due to transition to vibronic states at zero temperature, aparticular quantum state of the multi-mode molecular vibronic state maybe determined as an input state to the classical algorithm. Such aninput vibronic state may be a general M-mode vibronic state determinedand generated for an initial state of the molecule (e.g., a vacuum stateat zero temperature, a Gaussian state at a finite initial temperature,or Fock state, and the like). Such an input vibronic state correspondingto an initial vacuum state at zero temperature, for example, may be aGaussian state associated with the M vibronic modes. In some exampleimplementations, the input M-mode vibronic quantum state of the moleculeto the disclosed classical algorithm may be represented in a phase spaceas a quasi-probability distribution of a set of phase space variablesassociated with the phase space. Such a phase-space representation of aquantum state may be referred to as a P-representation.

An important property of a generalized P-representation may be that theunderlying phase-space distribution may always be chosen non-negativelysuch that the expectation value of a normal-ordered operator can bereadily estimated by sampling from the distribution. Such non-negativephase-space distribution may be referred to as positiveP-representation.

Unlike traditional classical algorithm for the calculation of molecularvibronic spectra, the quantum-inspired classical algorithm can providefast and efficient estimation of the molecular vibronic spectra. Thegain of such a classical algorithm over traditional classicalcalculation with complexity that grows combinatorially with respect tothe size of the problem or the number of the vibronic modes of themolecular system mostly relies on the principle that information beingsought or of interest (e.g., the molecular vibronic spectra) representsa coarse information that need not require detailed knowledge whosecalculation cannot be avoided when the traditional classical calculationmethod is being used. The example quantum-inspired algorithm describedbelow preforms sampling and analytics that more resemble the quantumsimulation in that it is not capable of revealing or calculating theseknowledge items pertinent to the molecular vibronic system, andfortunately, those knowledge items are not necessary and thus need notbe calculated for the estimation of the molecular vibronic spectra ofinterest.

As shown in 110 of FIG. 1 , this classical algorithm starts by preparingan input state of the M-mode molecular vibronic system based on aninitial vibronic state of interest (e.g., a vacuum state at zerotemperature, a Gaussian state at a finite temperature, a Fock state, andthe like). For an initial vacuum state, the input state may correspondto an Gaussian state, as described in further detail below. Merely as anexample, the input state may be expressed in the P-representation in aphase space. Again, the P-representation may be characterized by aquasi-probability distribution function of a bosonic state over thephase space that can always be chosen non-negative. For example thephase space variables for the quasi-probability distribution function ofthe P-representation may include a set of M-dimensional vectors. For aparticular example, the phase space variables may include a set of twoM-dimensional vectors. In FIG. 1 , such example quasi-probabilitydistribution function of the example input Gaussian state is shown by112 and represented by P(α, β), where α and β represent the set of twoM-dimensional vectors.

As further shown in 110 of FIG. 1 , the input state of the M-modemolecular vibronic state as represented by P(α, β) may then betransformed to generate a transformed P-representation P(α′,β′), asshown by 114, by applying an operator U, shown by 113, to the set of twophase space variables α and β to transform them into α′ and β′. Thephase space of α and β may be further sampled. Specifically, 2M complexnumbers of (α, β) may be sampled. The 2M samples may be transformed bythe operator U, as expressed by (α, β)→(α′, β′)=(Uα, Uβ). The operator Uis implemented as a matrix and represents the transmission/transitionmatrix for the M-mode molecular system in the phase space. The operatorU, for example, may comprise a mix of position and momentum operators.

As further shown by 115 of FIG. 1 , Fourier components of the molecularvibronic spectra may be estimated classically from the phase spacesampling above. Such estimate may be exact or may be approximated to acertain predetermined error level. For example, the Fourier componentsof the molecular vibronic spectra may by represented by {tilde over(G)}(k), where k represents the Fourier space of the molecular vibronicspectra. The classical process for estimating {tilde over (G)}(k) fromthe transformed sampled phase space variables (α′, β′) is provided infurther detail below for various initial and input states of the Mvibronic modes of the molecule. Such estimate may be based on acorrespondence between the Fourier components of the molecular vibronicspectra and an effective photon number operator {circumflex over (n)},as shown by 115 of FIG. 1 . The estimated Fourier components {tilde over(G)}(k) are exemplarily illustrated in 116 of FIG. 1 . Thereafter, theestimated Fourier components {tilde over (G)}(k) may then be used togenerate the molecular vibronic spectra 130 via a classical inverseFourier transform process, as indicated by 118 in FIG. 1 .

Although the example of FIG. 1 illustrates sampling in the Prepresentation phase space, the underlying principles apply toimplementations where the sampling is performed in other spaces. Otherrepresentations or spaces may be used under different circumstances,depending on, for example, the initial state of the molecular system.For example, the sampling space may be generally a mix of position andmomentum.

The example procedure above for estimating the molecular vibronicspectra in the quantum-inspired classical manner is further described inthe data and logic flow 200 of FIG. 2 . As shown in FIG. 2 , the dataand logic flow may include step 210 for determining a unitary operationconfigured to project an input superposition state of the bosonic systeminto number states corresponding to the M bosonic modes. Such an inputstate may be determined form an initial state of moledule system. Forexample, such an input state may be a Gaussian state of the M bosonicmodes corresponding to an vacuum initial state at zero temperature, forcalculating zero-temperature vibronic spectra of the molecule.

In step 220, each of N sampling bosonic states of the bosonic system maybe processed in a representation space of the bosonic system e.g., aphase space to generate N sets of samples based on at least the unitaryoperation, wherein each set of samples correspond to one of the Nsampling bosonic states and are sampled from at least one variable inthe representation space, N being a positive integer.

In step 230, Fourier components the transition spectra may beanalytically generated based on the N sets of samples. Such generationmay be exact or may be approximate to a predetermined error levelcorresponding to the size (number of vibronic modes) of the molecule.

In step 240, inverse-transformation may be performed on the Fouriercomponents to estimate the transition spectra. The correspondencebetween the Fourier components and the final transition spectra may beexact or may be approximate to a predetermined error level correspondingto the size (number of vibronic modes) of the molecule for differentinitial and input sates, as described in further detail below.

The example implementations above may be applied to generate vibronicspectral for a molecular system at zero temperature. In particular, theinitial vibronic state of the molecular system would be the vacuum stateof all M modes. The corresponding input state to the sampling processabove may be an Gaussian state of the M vibronic modes, and thecorrespondence between the Fourier components and the final transitionspectra may be exact, as described in further detail below. Therepresentation space in such a case may be the P-representation space.

In some other specific example implementations, vibronic spectra may beestimated for a molecular system at a finite non-zero temperaturefollowing the zero-temperature implementation. For example, as describedin further detail below, M auxiliary modes may be added in comparison tothe zero-temperature case. Two-mode squeezed states may be prepared tosimulate thermal states and detect the additional auxiliary modes aswell, as depicted in FIG. 4 and described in further detail below. Assuch, the procedure described above and in the data and logic flow ofFIG. 2 for estimating zero-temperature vibronic spectra of the molecularsystem applies to the finite temperature situation, except that, as acost, the additional M auxiliary modes are included in the preparationof the input Gaussian state and in the sampling in the phase space.

In yet some other example implementations, the sampling method above maybe applied to an estimation of molecular vibronic spectra fornon-Gaussian initial state, such as a Fock state. In order to employ amapping from molecular vibronic spectra to Gaussian boson sampling,additional M modes and two-mode squeezed states of some squeezingparameters may be introduced. In comparison to the finite temperatureGaussian initial state case, a post selection may be further performedfor the additional modes. As described in further detail below, theestimation of the Fourier components may not be exact in such asituation. The error may exponentially increase with the photon numberin the initial Fock state, suggesting an exponential increase in samplesneeded according to the sampling process of FIGS. 1 and 2 . However, asdescribed in further detail below, the above sampling process may becombined with, for example, the Gurvits' algorithm generalized toinclude multiphotons to approximate the Fourier components of themolecular spectra with reduced costs in such situation within areasonable error. The sampling may be performed in the particularlydetermined representation space as described in further detail below.

The present disclosure implies that Gaussian boson sampling itself is ahard task, but when it comes to the coarse-grained version of sampling,an efficient classical counterpart may exist. Therefore, on the contraryto the common understanding, the disclosure herein suggests thatmolecular vibronic spectra may not be the candidate for which only thepower of a quantum sampler can boost the computational performancebeyond classical means. Rather, classical algorithm as described abovemay diminish or reduce the quantum advantage.

On the other hand and as shown in further detail below, beyond Gaussianframework where a direct integral does not work, the disclosure hereinshows that a Monte Carlo method may be inefficient when the number ofquanta of the input state is large. When a similar problem insingle-photon boson sampling is considered, computing the Fouriercomponents is a #P-hard problem, which is in stark contrast to Gaussianboson sampling case. Nevertheless, a quantum-inspired classicalefficient algorithm has been identified herein to approximate Fouriercomponents of molecular vibronic spectra, which may be further appliedto other application scenarios.

Non-Sampling Classical Molecular Vibronic Spectra Computation

The transition probability between an initial vibronic mode and acertain final vibronic mode is referred to as the Franck-Condon (FC)factor. FC profiles (FCP) may be obtained by computing many FC factorscorresponding to a given vibrational transition frequency (ω_(vib)). Forexample, the FC factor at 0 K may be obtained with the initial vacuumstate as

$\begin{matrix}{{{{FCP}\left( \omega_{vib} \right)} = {\sum\limits_{m = 0}^{\infty}{{❘\left\langle {m{❘{\hat{U}}_{Dok}❘}0} \right\rangle ❘}^{2}{\delta\left( {\omega_{vib} - {\omega^{\prime} \cdot m}} \right)}}}},} & (1)\end{matrix}$

where the Doktorov transformation Û_(Dok) is given by

Û _(Dok) ={circumflex over (D)}(δ/√{square root over (2)})S _(Ω) ^(\),{circumflex over (R)} _(U) Ŝ _(Ω),  (2)

where δ(·) represents the delta function and m=(m₁, . . . , m_(M)) isthe final vibronic modes excitation vector with M vibronic modes. Here,Ŝ, {circumflex over (D)}, and {circumflex over (R)}_(u) representsqueezing, displacement, and rotation operators, respectively, and Ω andΩ′ are given by

Ω≡diag(√{square root over (ω₁)}, . . . ,√{square root over(ω_(M))}),Ω′≡diag(√{square root over (ω′₁)}, . . . ,√{square root over(ω′_(M))}),  (3)

where ω_(i)'s and ω_(i)'s account for initial and final harmonic angularfrequencies, respectively.

The FCP above represents the molecular vibronic spectra. Eq. (1) showsthat the direct computation of all FC factors are inevitablycomputational hard. In particular, not only is the number ofcombinations of m satisfying the delta function combinatorially large,the FC factor |

m|Û_(Dok)|0

|² corresponds to the hafnian of a matrix, the computation of whichcosts exponentially in the size of matrix or a number of vibronic modes.For a fixed m=Σ_(i=1) ^(M)m_(i), the number of FC factors may be givenby

$\begin{matrix}{\begin{pmatrix}{M + m - 1} \\m\end{pmatrix}.} & (4)\end{matrix}$

Gaussian Boson Sampling for Molecular Vibronic Spectra Generation viaQuantum Simulation

Gaussian boson sampling represents one of the sampling problems that areproved hard for classical computers under some plausible assumptions. Toimplement Gaussian boson sampling, a product of M single-mode Gaussianstates, |ψ_(in)

, may first be prepared, and the product state may be fed into an M-modelinear-optical circuit Û_(LON), composed of beam splitters. The outputstate |ψ_(out)=Û_(LON)|ψ_(in)

may be detected or measured by photon-number detectors.

The output photon patterns in may follow the probability distributionbelow

ρ(m)=|

m|Û _(LON)|ψ_(in)

|²,  (5)

which may be related to the (loop) hafnian of a certain matrix.

Due to Gaussian boson samplings computational power beyond classicalcomputers, a Gaussian boson sampler may be employed to generatemolecular vibronic spectra. For example, the correspondence between theFC factor |

m|Û_(Dok)|0

|² in Eq. (1) for zero temperature and the output probability ρ(m) ofdetection photon number pattern in in Gaussian boson sampling may beestablished with an appropriate choice of input state and thelinear-optical operation.

In some implementations, an appropriate Gaussian boson sampling mayalways be chosen by exploiting the fact that the Doktorov operationÛ_(Dok) is a Gaussian unitary, and can be decomposed asÛ_(LON){circumflex over (D)}(α)Ŝ(r){circumflex over (V)}_(LON). Since{circumflex over (V)}_(LON)|0

=|0

, the decomposed unitary operation can be implemented by preparingdisplaced squeezed states and applying a linear-optical circuit, whichis exactly the procedure of Gaussian boson sampling. In other words, forthe Doktorov operation Û_(Dok) of a molecular structure, a Û_(LON) canbe identified and used to build the linear optical circuits. The linearoperators above, may be mix of both position and momentum operators.

Therefore, by implementing Gaussian boson sampling and post-processingthe output samples using the relation (1), FCP(ω_(vib)) may beefficiently generated via simulation by the linear optical quantumcircuit, as shown by 120 of FIG. 1 and as described above.

More specifically, sampling may be first taken from ρ(m) using Gaussianboson sampler. A number of samples may be collected, represented byS={m^((i))}_(i=1) ^(N). The samples may then be sorted according to thedelta function above with respect to the frequency (or energy). In someimplementations, a window around each frequency Δω_(vib) may beintroduced. The window may be provided with a predetermined and/orconfigurable size Δω_(vib).

In other words, when a samples frequency ω′·m lies on [ω_(vib)Δω_(vib)/2,ω_(vib)+Δω_(vib)/2] for a certain ω_(vib), a statistical bincorresponding to ω_(vib) may be incremented. Once sufficient number ofsamples are post-processed, an estimate of eFCP(ω_(vib)) (denotingempirical FCP, referring to the experimental quantum simulation aspectof the methodology) may be generated.

Exact Quantum-Inspired Classical Solution for Molecular Vibronic Spectraat Zero Temperature

Inspired by the mapping from a molecular vibronic spectra problem to aGaussian boson sampling and the phase-space method in quantum optics,the quantum-inspired solution for classically obtaining the molecularvibronic spectral without replying on a quantum mechanical simulationmachine may be identified.

For example, it is noted that the Fourier components of the molecularvibronic spectra may be analytically derived and simplified as:

$\begin{matrix}{{\overset{\sim}{G}(k)} \equiv {\left\langle {\underset{i = 1}{\overset{M}{\otimes}}{\sum\limits_{m_{i} = 0}^{\infty}{❘m_{i}}}} \right\rangle\left\langle {m_{i}❘_{i}e^{{- {ikm}_{i}}\omega_{i}^{\prime}\theta}} \right\rangle}} & (6)\end{matrix}$ $\begin{matrix}{= \left\langle e^{{- {ik}}{\hat{n} \cdot \omega_{i}^{\prime}}\theta} \right\rangle} & (7)\end{matrix}$ $\begin{matrix}{{G\left( \omega_{vib} \right)} = {\frac{1}{K + 1}{\sum\limits_{k}{{\overset{\sim}{G}(k)}e^{{ik}\theta\omega_{vib}}}}}} & (8)\end{matrix}$ $\begin{matrix}{{= {\sum\limits_{m = 0}^{\infty}{{p(m)}{\delta\left( {\omega_{vib} - {\omega^{\prime} \cdot m}} \right)}}}},} & (9)\end{matrix}$

where θ≡2π/(K+1), k∈{0, . . . , K}, K is the number of bins indiscretized distribution of ω_(vib), and |m_(i)

m_(i)|_(i), is the projector on m_(i), photon state for the ith mode.

The final expression above shows the exact correspondence betweenG(ω_(vib)) and FCP (ω_(vib)). Hence, in general, a Fourier component isan overlap between an output state |ψ_(out)

and a phase-shifted state e^(−ikñ·w′θ)|ψ_(out)

. Once such overlap is computed efficiently (e.g., in a classicalmanner), the Fourier components may then be obtained and used to recoverthe spectra by an inverse Fourier transform.

Especially when the relevant quantum states are Gaussian states, whichis the case for molecular vibronic spectra at zero temperature totransition to thermal states, the Fourier components can be analyticallycomputed. One example way to compute the Fourier components may involveusing the positive P-representation of the Gaussian state in a phasespace as described above.

More generally, a generalized P-representation of an M-mode quantumstate may represent one of the quasi-probability distributions of acorresponding bosonic state of the molecular vibronic modes:

$\begin{matrix}{{\hat{\rho} = {\int_{{\mathbb{C}}^{2M}}{{P\left( {\alpha,\beta} \right)}{\hat{\Lambda}\left( {\alpha,\beta} \right)}d\alpha d\beta}}},{{\hat{\Lambda}\left( {\alpha,\beta} \right)} \equiv \frac{{\left. {❘\alpha} \right\rangle\left\langle \beta^{*} \right.}❘}{\left\langle \beta^{*} \middle| \alpha \right\rangle}},} & (10)\end{matrix}$

where |α

=|α₁

⊗ . . . ⊗|α_(M)

is an M-mode coherent state.

An important property of the generalized P-representation is that thedistribution P(α, β) representing the bosonic state can always be chosennon-negative (hence positive P-representation), and the expectationvalue of a normal-ordered operator can be readily computed.Specifically, for molecular vibronic spectra generation, 2M complexnumbers (α,β) may be sampled from the generalized P-representation of aninput Gaussian state and the samples may be transformed by thelinear-optical unitary operation U as (α,β)→(α′,β′)=(Uα, U*β), where Urepresents the trans-mission/transition matrix in the phase space. Forthe measurement portion of the analytics, the photon number operator{circumflex over (n)}_(i) may be replaced with α′β′. For example, forthe measurement outcome (m₁, m₂) for two modes, the correspondingobservable may be written as |m₁

m₁|₁⊗|m₂

m₂|₂=: {circumflex over (n)}₁ ^(m) ¹ e^(−{circumflex over (n)}) ¹/m₁!⊗{circumflex over (n)}₂ ^(m) ² e^(−{circumflex over (n)}) ² /m₂!:,so that the phase-space variable represented by (α′₁β′₁)^(m) ¹(α′₂β′₂)^(m) ² e^(−α′) ¹ ^(β′) ¹ ^(−α′) ² ^(β′) ² /m₁!m₂! may beobtained. Here, |j

j|_(i) represents the projector on j photon state for the ith mode, and:Ô: represents the normally ordered form of an operator Ô. For asufficiently large number of samples, it may be guaranteed that fornormal-ordered observables, the average converges to the expectationvalue for the output quantum state

For example, the positive P-representation for a single-mode squeezedstate may be chosen as:

$\begin{matrix}{{{P\left( {\alpha,\beta} \right)} = {\frac{\sqrt{1 + \gamma}}{\pi\gamma}e^{{{- {({\alpha^{2} + \beta^{2}})}}{({\gamma^{- 1} + {1/2}})}} + {\alpha\beta}}}},} & (11)\end{matrix}$

where the domain is the real plane and γ≡e^(2r)−1 for a squeezingparameter r>0.

Using the positive P-representation of single-mode Gaussian states andthe relation for a normal-ordered operator such ase^(−γ{circumflex over (n)})=: e^({circumflex over (n)}(e) ^(−γ) ⁻¹⁾:,the Fourier components may be rewritten as:

$\begin{matrix}{{\overset{\sim}{G}(k)} = \left\langle {:{\exp\left\lbrack {\sum\limits_{i = 1}^{M}{{\hat{n}}_{i}\left( {e^{{- {ik}}\omega_{i}^{\prime}\theta} - 1} \right)}} \right\rbrack}:} \right\rangle} & (12)\end{matrix}$ $\begin{matrix}{{= {\int_{{\mathbb{R}}^{2M}}{d\alpha d\beta{P_{in}\left( {\alpha,\beta} \right)}{\exp\left\lbrack {\sum\limits_{i = 1}^{M}{\alpha_{i}^{\prime}{\beta_{i}^{\prime}\left( {e^{{- {ik}}\omega_{i}^{\prime}\theta} - 1} \right)}}} \right\rbrack}}}},} & (13)\end{matrix}$

where P_(in)(α,β) is the positive P-representation of an input squeezedstate. Here, :Ô: represents the normal-ordered form of an operator Ô,and (α′,β′)=(Uα,U*β) accounts for the linear-optical unitary operation.

More specifically, when a coherent state is processed through an M-modelinear-optical circuit, it is transformed as Û|α

=|Uα

, where U is the corresponding M×M unitary matrix for the circuit.

Since {acute over (G)}(k) is now written as a Gaussian integral, it canbe analytically obtained, which is given by

$\begin{matrix}{{{\overset{\sim}{G}(k)} = {\mathcal{N}\frac{\left( {2\pi} \right)^{M}}{\sqrt{\det(Q)}}{\exp\left( {{\frac{1}{2}c^{T}Q^{- 1}c} + c_{0}} \right)}}},} & (14)\end{matrix}$ whereϕ_(j) = −kθω_(j), Q ≡ ( 2 ⁢ Γ - 1 + M - U T ⁢ diag ⁡ ( ei ⁢ ϕ j ) j = 1 M ⁢ U * - U † ⁢ diag ⁡ ( e i ⁢ ϕ j ) j = 1 M ⁢ U 2 ⁢ Γ - 1 + M) , and , ( 15 ) $\begin{matrix}{{\mathcal{N} \equiv {\prod\limits_{i = 1}^{M}\frac{\sqrt{1 + \gamma_{i}}}{\pi\gamma_{i}}}},{\Gamma \equiv {{diag}\left( \gamma_{i} \right)}_{i = 1}^{M}},{\Phi = {{diag}\left( {e^{i\phi_{j}} - 1} \right)}_{j = 1}^{M}},} & (16)\end{matrix}$${c \equiv \left( {a^{T},b^{T}} \right)^{T}},{a \equiv \frac{U^{T}\Phi\delta}{\sqrt{2}}},{b \equiv \frac{U^{T}\Phi\delta}{\sqrt{2}}},{c_{0} \equiv {\frac{\delta^{T}\Phi\delta}{2}.}}$

Here, δ represents a real vector accounting for the displacement in theDoktorov transformation (see Eq. (2)). More detailed derivation aboveand proof of the convergence of the integral are further provided laterin this disclosure.

It implies that Fourier components of the molecular vibronic spectra atzero temperature have a solution. As such, the spectra may be obtainedby simply taking the inverse Fourier transform, thereby avoiding havingto use a quantum simulator to generate the molecular vibronic spectra.

As an example, FIG. 3 shows a comparison between an analytical solutionof an example molecular vibronic spectra based on the formalism aboveand direct classical computation. FIG. 3 shows that the analyticalsolution based on the Fourier formalism above matches with the spectraobtained by directly computing probabilities. For the analyticalsolution, the Fourier components are first computed using Eq. (14)followed by the inverse Fourier transform to obtain the spectra. Thecorresponding Doktorov transformation is determined in accordance withthe simulated molecular system being CH₂O₂, 1¹A′→1²A′. The number ofsamples processed in FIG. 2 is N=10⁴. An increase of number of bins inthe frequency domain may be achieved with a polynomial computationalcost.

Exact Solution for Molecular Vibronic Spectra at Finite Temperature

The same method above may be applied to initial states of the molecularsystem at finite temperatures. Specifically, in some exampleimplementations, the quantum-inspired classical approach to thegeneration of molecular vibronic spectra may be further applied tofinite temperatures other than zero temperatures.

At finite temperature T, the molecular vibronic spectra may be writtenas

$\begin{matrix}{{{{FCP}\left( \omega_{vib} \right)} = {\sum\limits_{n,{m = 0}}^{\infty}{{p\left( {n,m} \right)}{\delta\left( {\omega_{vib} - \left( {{\omega^{\prime} \cdot m} - {\omega \cdot n}} \right)} \right)}}}},{where}} & (17)\end{matrix}$ $\begin{matrix}{{p\left( {n,m} \right)} = {{p_{T}(n)}{{❘\left\langle {n{❘{\hat{U}}_{Dok}❘}m} \right\rangle ❘}^{2}.}}} & (18)\end{matrix}$

Here, p_(T)(n)≡

n|{circumflex over (ρ)}_(T)|n

is the probability of photon number |n

from an M-mode thermal state {circumflex over (ρ)}_(T) at temperature T.

In some specific implementations, a Gaussian boson sampling quantumprocedure described above may be correspondingly constructed for thefinite temperature situations. For example, M auxiliary modes may bedetermined and added to the quantum simulation. Two-mode squeezed statesmay be further prepared to simulate thermal states. The photon numbersin the additional modes may be detect as well. The process is exemplarydepicted in FIG. 4 . Inspired by such quantum simulation, a classicalsampling algorithm of Fourier components in the P-representation similarto the zero-temperature situation may also be employed with an additionof the auxiliary modes.

In such example implementations, the original M modes output photonscorrespond to in and the additional M modes' output photons correspondsto n. Therefore, the same framework described above for zero temperaturemay be employed to obtain the solution for molecular vibronic spectra atfinite temperature, as illustrated in FIG. 4 , showing a comparison ofstate transform between zero-temperature scenario and the finitetemperature scenario.

For example, the FCP of the molecular system at the finite temperaturemay be rewritten as:

$\begin{matrix}{{{FCP} = {\sum\limits_{m,{n = 0}}^{\infty}{{❘\left\langle {m,{n{❘{\hat{U}}_{Dok}^{\prime}❘}0}} \right\rangle ❘}^{2}{\delta\left( {\omega_{vib} - \left( {{\omega^{\prime} \cdot m} - {\omega \cdot n}} \right)} \right)}}}},} & (19)\end{matrix}$

which is the same form as Eq. (1) except for the negative part in thedelta function. Here, Û′_(Dok)=Û_(Dok)Û_(TMSO) is obtained by addingtwo-mode squeezing operators. The term “TMSO” denotes“Two-Mode-Squeezing Operator”. Each pair of two modes include anoriginal mode and an auxiliary mode.

Again, the Fourier transform of the FCP at finite temperature may begiven by

{tilde over (G)}(k)=

e ^(−ik{circumflex over (n)}·ω″θ)

,  (20)

where ω″=(ω′^(T), −ω^(T))^(T). Therefore, the analytic solution forFourier components can thus be found for the molecular system at afinite temperature. The vibronic spectra at the finite temperature thusmay be estimated by inverse Fourier transform of the analyticallygenerated Fourier components above.

Molecular Vibronic Spectra for Non-Gaussian Initial State

The above implementations assume a vacuum or otherwise a thermal orGaussian initial state for the molecular vibronic modes. In a moregeneral case, for example, the molecular system may be in a non-thermalinitial state. In such circumstances, the molecular vibronic system maynot be described by Gaussian states.

In a particular circumstance, the initial state of the molecularvibronic system may be a Fock state or number state of the vibronic, |n

, rather than vacuum |0

a Gaussian state. For optical processes including the single vibroniclevel fluorescence and the resonance Raman scattering, for example, aspectra from specific vibronic levels needs to be analyzed, instead of athermal distribution or a ground state. Indeed, a quantum simulation ofsuch a process has also been experimentally implemented for variousphotoelectron processes such as photodetachment of ozone anion. In thesecases, the FCP of the molecular vibronic system for such a non-Gaussianinitial state may be expressed as:

$\begin{matrix}{{{{FCP}\left( \omega_{vib} \right)} = {\sum\limits_{m = 0}^{\infty}{{❘\left\langle {m{❘{\hat{U}}_{Dok}❘}n} \right\rangle ❘}^{2}{\delta\left( {\omega_{vib} - {\omega^{\prime} \cdot m}} \right)}}}},} & (21)\end{matrix}$ $\begin{matrix}{{{\overset{\sim}{G}(k)} = {\left\langle {n{❘{{\hat{U}}_{Dok}^{\dagger}e^{{- {ik}}{\hat{n} \cdot \omega^{\prime}}\theta}{\hat{U}}_{Dok}}❘}n} \right\rangle = \left\langle {n{❘\hat{V}❘}n} \right\rangle}},} & (22)\end{matrix}$

where {circumflex over (V)}≡Û_(Dok)^(†)e^(−k{circumflex over (n)}·w′θ)Û_(Dok) is a Gaussian unitaryoperation. Using Bloch-Messiah decomposition, the Gaussian unitaryoperation is decomposed as {circumflex over (V)}={circumflex over(D)}(ξ)Û_(lin2)Ŝ(r)Û_(lin1). Because of the non-Gaussian initial state,the spectra cannot be directly mapped to a Gaussian boson sampling andso it may require a non-Gaussian type of boson sampling.

In some implementations, the Fourier components can be written as a loophafnian of a matrix:

$\begin{matrix}{{{\overset{\sim}{G}(k)} = \frac{{lhaf}\left( {\overset{\sim}{\Sigma}}_{n} \right)}{{n!}Z}},{where}} & (23)\end{matrix}$ $\begin{matrix}{\Sigma \equiv \begin{pmatrix}{U_{{lin}2}\tanh{rU}_{{lin}2}^{T}} & {U_{{lin}2}{sech}{rU}_{{lin}1}} \\{U_{{lin}1}^{T}{sech}{rU}_{{lin}2}^{T}} & {{- U_{{lin}1}^{T}}\tanh{rU}_{{lin}1}}\end{pmatrix}} & (24)\end{matrix}$

is an 2M×2M complex symmetric matrix,

$\begin{matrix}{Ϛ \equiv \begin{pmatrix}{{{- U_{{lin}2}}\tanh{{rU}_{{lin}2}^{T} \cdot \xi^{*}}} + \xi} \\{{- U_{{lin}1}^{T}}{sech}r{U_{{lin}2}^{T} \cdot \xi^{*}}}\end{pmatrix}} & (25)\end{matrix}$

is an 2M-dimensional vector, and

$\begin{matrix}{{{Z^{- 1} \equiv \left\langle {0{❘\hat{V}❘}0} \right\rangle} = \frac{e^{\frac{1}{2}{\xi^{T} \cdot \Sigma \cdot \xi}}}{\sqrt{\prod_{i = 1}^{M}{\cosh r_{i}}}}},} & (26)\end{matrix}$

is the normalization factor, which is the same as the Fourier componentsof molecular vibronic spectra at zero temperature.

Here, {tilde over (Σ)}_(n) is obtained by first replacing the diagonalelements of Σ by ζ to obtain {tilde over (Σ)} and repeating ith row andcolumn of each block matrix of {circumflex over (Σ)}n_(i), times, and isthus an n×n matrix with n≡Σ_(i=1) ^(M)n_(i). Therefore, computing theFourier components reduces to computing loop hafnians of n×n complexsymmetric matrices.

Loop hafnian is a quantity that is related to perfect matchings of agraph including loops (hafnian does not allow loops). The best-knownalgorithms computational cost of computing a loop hafnian isO(n³2^(n/2)) with n being the matrix size. Thus, if FC factors, which iswritten as a loop hafnian of a matrix whose size is n+m is to bedirectly computed to obtain FC profiles in Eq. (1), then the complexitywould be exponential in (n+m)/2. On the other hand, the complexity ofcomputing the Fourier components in the proposed method relies on theinput photons n rather than the output photons m and the system size,i.e., the proposed method can be efficiently implemented as long as n issmall enough. More precisely, since the redundancy of rows and columnsfor n_(i)≥2 does not increase the complexity of computing a loop hafnianbecause it does not increase the rank of the matrix, the importantfactor for complexity is the number of nonzero elements of n.

Meanwhile, computing the loop hafnian of a general complex symmetricmatrix in multiplicative error is known to be #P-hard. It implies thatcomputing the Fourier components of general molecular vibronic spectralwith Fock state inputs and squeezing is also #P-hard. More formally, anarbitrary complex matrix may be embedded into a submatrix of a unitarymatrix with an appropriate normalization so the corresponding spectrageneration problem reduces to computing a permanent, which is #P-hard.The unitary matrix for Fourier coefficient is written asÛe^(i{circumflex over (n)}·ω)Û^(†), which is a diagonalized form ofarbitrary unitary matrix, as shown further below.

Therefore, the Fourier components for molecular vibronic spectra mappedto Gaussian boson sampling have an analytic solution while thecorresponding problem for non-Gaussian cases is a computationally hardproblem, which presents a separation between Gaussian boson sampling andnon-Gaussian type of boson sampling such as single-photon bosonsampling. The Fourier components reduce to hafnian and permanent whenthere is no displacement or squeezing, respectively.

Even if a boson sampling experiment is performed using a quantum devicethe spectra is reproduced from the outcomes, estimation of the spectrahas an additive error depending on the number of samples. Therefore,reproducing the spectra within a multiplicative error is not expected tobe achieved by running a boson sampling circuit, so an additive error isthe relevant quantity using a classical simulation to compare with aquantum device.

Fock States: Sampling Using Positive P-Representation

To employ a mapping from molecular vibronic spectra to Gaussian bosonsampling, again, additional M modes and two-mode squeezed states ofsqueezing parameters si may be introduced, as illustrated in FIG. 5 ,which shows an comparison between Gaussian boson sampling for a thermalstate, and mapping to Gaussian sampling for a Fock state. The differencefrom the finite temperature case above is that now the outcome n for theadditional modes need to be post-selected.

Thus, the corresponding FCP can be written as

$\begin{matrix}{{{G\left( \omega_{vib} \right)} = {\frac{1}{N(n)}{\sum\limits_{m = 0}^{\infty}{{p\left( {m,n} \right)}{\delta\left( {\omega_{vib} - {\omega^{\prime} \cdot m}} \right)}}}}},} & (27)\end{matrix}$

where

is the normalization due to the post-selection,

$\begin{matrix}{{{\mathcal{N}(n)} = {\prod\limits_{i = 1}^{M}\frac{\tanh^{2n_{i}}s_{i}}{\cosh^{2}s_{i}}}},} & (28)\end{matrix}$

and p(m, n)=

(n)

m|Û_(Dok)|n

².

Its Fourier components are given by

$\begin{matrix}{{\overset{\sim}{G}(k)} = {\frac{1}{\mathcal{N}}\left\langle {e^{{- {ik}}\theta{\sum_{j = 1}^{M}{{\hat{n}}_{j}\omega_{j}^{\prime}}}} \otimes {❘n}} \right\rangle{\left\langle \left. n \right| \right\rangle.}}} & (29)\end{matrix}$

Approximation of Hafnian

In some implementations, another method is to use the expression in Eq.(23). Assume that the displacement is zero, the loop hafnian in Eq. (23)then reduces to a hafnian. In this case, a randomized algorithm isconstructed to estimate a hafnian of an n×n complex symmetric matrix Σusing the equality

$\begin{matrix}{{{{haf}(\Sigma)} = {{\mathbb{E}}_{\upsilon \in {\{{0,1}\}}^{n}}\left\lbrack {\left( {- 1} \right)^{\sum\limits_{i = 1}^{n}\upsilon_{i}}\frac{2^{n/2}}{\left( {n/2} \right)!}\left( {h^{T}\Sigma h} \right)^{n/2}} \right\rbrack}},} & (30)\end{matrix}$

where

_(Σ)[·] is the average over the normal distribution of covariance matrixΣ, and h=(½−v₁, . . . , ½−v_(M))^(T) and n is even (otherwise thehafnian is zero).

Therefore, by sampling v∈{0,1}^(n) uniformly and averaging over thesamples, the hafnian may be estimated, which, together with the factthat ∥Σ∥=1 in Eq. (24), leads to the bound for estimate q of a Fouriercomponent as

$\begin{matrix}{{{\Pr\left\lbrack {{❘{q - \frac{{haf}\left( \Sigma_{n} \right)}{Z{n!}}}❘} \gtrsim {\epsilon\frac{e^{n/2}}{\sqrt{\pi n}Z}}} \right\rbrack} \leq {2e^{{- N}{\epsilon^{2}/2}}}},} & (31)\end{matrix}$

where Z=√{square root over (Π_(i=1) ^(n) cosh r_(i))} and N is thenumber of samples. Thus, if Π_(i=1) ^(n) cosh r_(i)>e^(n), theestimation error is smaller than e with exponentially small failureprobability if the number of samples is chosen as N=O(1/ϵ²).

It suggests that if the condition is satisfied, classically estimatingthe Fourier components and taking the inverse Fourier transformationrenders the same scaling of precision to running a boson samplingcircuit.

For nonzero displacement and loop hafnian, although there is a similarequality as in Eq. (30), the resultant error bound is more complicatedthan hafnian.

It is further provided another method of estimating a loop hafnian thatuses the integral representation of loop hafnian in Methods. Since loophafnian corresponds to the high order moment of multivariate normaldistribution, loop hafnian of an n×n complex symmetric matrix {tildeover (Σ)} may be written as an integral:

$\begin{matrix}{{{lhaf}\left( \overset{\sim}{\Sigma} \right)} = {\frac{1}{\pi^{4n}}{\int_{{\mathbb{C}}^{4n}}{d\alpha d\gamma d\omega d\beta{\underset{i = 1}{\overset{n}{\prod}}{{\left\lbrack {\left( {\alpha_{i} + \xi_{i} + \zeta_{i}^{\alpha}} \right){\left( {\beta_{i} + \zeta_{i}^{\beta^{*}}} \right)^{*}}} \right\rbrack \times {{\exp\left( {{- \frac{1}{2}}{{\upsilon\left( {\alpha,\beta,\gamma,\omega} \right)}^{T} \cdot \Sigma^{- 1} \cdot {\upsilon\left( {\alpha,\beta,\gamma,\omega} \right)}}} \right)}.}}}}}}}} & (32)\end{matrix}$

The expression is obtained in the section below entitled “Approximationalgorithm of hafnian’,’ where the explicit definition of each term isgiven. Therefore, the loop hafnian can be estimated by decomposing theexponential term into an imaginary part and a real part, which followsthe normal distribution, and, by sampling from the normal distributionand averaging over the random variable including the imaginary part.

Fock States: Approximation of Fourier Components of Molecular VibronicSpectra

The observation that the additive error is the achievable error fromrunning a quantum device, some classical algorithms that might bepotentially able to efficiently estimate the Fourier components withinan additive error may be conceived.

For example, Gurvits's algorithm can efficiently approximate thecorresponding quantity within an additive error e in single-photon bosonsampling, where the FCP is written as permanent of a submatrix of aunitary matrix (see the section entitled “Computing Fourier componentsof vibronic spectra for single-photon input cases” below). Indeed, if aboson sampling is directly run to obtain the FCP, the Chernoff boundposes a limitation of accuracy, which is given by O(1/poly(N)) when thenumber of samples is N. It suggests that even though exactly computingthe Fourier components might be difficult, it may be possible toefficiently approximate them within a reasonable additive error.

Since inverse Fourier transformation does not amplify the additive errorof Fourier components, if there is an efficient algorithm thatapproximates the Fourier components within an error e in a running timeO(poly(1/ϵ)), the algorithm may allow for achieving a similar accuracyas running a boson sampling.

Fock-State: Estimating Fourier Components of Spectra with Boson Samplingwith Generalized Gurvits' Algorithm

For Fock-state boson sampling, the Fourier components can be written as

$\begin{matrix}{{{\overset{\sim}{G}(k)} = {{\left\langle {n{❘{{\hat{U}}^{\dagger}e^{{- {ik}}\theta{\hat{n} \cdot \omega}}\hat{U}}❘}n} \right\rangle \equiv \left\langle {n{❘V❘}n} \right\rangle} = \frac{{Per}\left( V_{n,n} \right)}{n!}}},} & (33)\end{matrix}$

where a unitary operator {circumflex over (V)} is defined as {circumflexover (V)}≡Û^(†)e^(−k{circumflex over (n)}·ωθ)Û, consisting of Û and aphase-shift operator e^(−ik{circumflex over (n)}·wθ) and Û^(†).

Thus, it is a linear-optical circuit characterized by a unitary matrixV=U^(†)DU, with D≡diag(e^(−ikθω) ¹ , . . . , e^(−ikθω) ^(M) )characterizing the phase-shift operator. Here, such a diagonal formsuggests that V is a general form of an arbitrary unitary matrix.Together with this fact, since computing the permanent of a generalmatrix in multiplicative error is #P-hard and one can embed an arbitrarycomplex matrix into a submatrix of a unitary matrix by normalizing thematrix, computing its Fourier component in multiplicative error is also#P-hard. Therefore, it shows that computing Fourier components inmultiplicative error is a #P-hard problem, and consequently, it provesthat the exact computation of the spectra is also a #P-hard problem.

In fact, it is noted that even if a boson sampling experiment runs usinga quantum device and reproduces the spectra from the sampling outcomes,the resultant spectra has an additive sampling error. Therefore,reproducing the spectra within a multiplicative error is not expected tobe achievable by running a boson sampling circuit with polynomially manysamples. As such, an additive error is a relevant target using aclassical simulation to compare with a quantum device.

Interestingly, there exists a classical algorithm, the so-calledGurvits' algorithm, which can efficiently approximate the permanent of amatrix within an additive error. However, this algorithm is notsufficient for general Fock-state boson sampling cases where n_(i)contains more than a single-photon for some i's because the error canincrease exponentially in the system size in that case (described infurther detail below). The Gurvits' algorithm, however, may be modifiedto estimate the spectra which includes multiphotons n using thefollowing equality (see detailed derivation in the next section):

$\begin{matrix}{{\frac{{Per}\left( V_{n,n} \right)}{n!} = {{\mathbb{E}}_{x \in \mathcal{X}}\left\lbrack {\prod\limits_{i = 1}^{M}\left( \frac{{{\overset{\_}{y}}_{i}({Vy})}_{i}}{n_{i}} \right)^{n_{i}}} \right\rbrack}},} & (34)\end{matrix}$

where y_(i)≡√{square root over (n_(i))}x_(i), x∈

≡

[n₁+1]× . . . ×

[n_(M)+1], where

[j] is the set of jth roots of unity.

Thus, by sampling the random variable with uniform x∈

, the randomized algorithm gives an estimate μ of permanent of an n×nmatrix such that

$\begin{matrix}{{{❘{\frac{{Per}\left( V_{n,n} \right)}{n!} - \mu}❘} < {\epsilon{V}^{n}}},} & (35)\end{matrix}$

with high probability 1−δ with running time T=O(poly(n, 1/ϵ, log δ⁻¹)).Here, ∥V∥ is the spectral norm of the matrix V, which is always 1 forthis case because V is a unitary matrix.

Thus, the modified (or generalized) Gurvits' algorithm enables anefficiently approximation of Fourier components within a reasonableadditive error even though computation of the Fourier components in amultiplicative error is hard (#P-hard). Note that although a modified orgeneralized Gurvits' algorithm might be used to approximate individualprobabilities in Eq. (9), it is nontrivial to approximate the sum ofexponentially many probabilities, which is enabled by approximatingFourier components instead of probabilities.

The remaining challenge is the propagation of the error of Fouriercoefficients to that of the spectra through the inverse Fouriertransformation. Using Parseval's relation, we prove that as long as weestimate the Fourier coefficients with a small error ϵ, the transformedspectra's error is also small as ϵ:

$\begin{matrix}{{{\sum\limits_{\Omega = 0}^{\Omega_{\max}}{❘{\Delta{G(\Omega)}}❘}^{2}} = {{\frac{1}{\Omega_{\max} + 1}{\sum\limits_{k = 0}^{\Omega_{\max}}{❘{\Delta{G(k)}}❘}^{2}}} \leq \epsilon^{2}}},} & (36)\end{matrix}$

which proves that |ΔG(Ω)|≤ϵ for any Ω, where ΔG and Δ{tilde over (G)}represent the error of spectra and Fourier component estimation.

Hence, if there is an efficient algorithm that approximates the Fouriercomponents within an error e in a running time T=O(poly(M, 1/ϵ, logδ⁻¹)), the algorithm enables us to achieve the same accuracy as runninga boson sampling. For Fock-state boson sampling case, Gurvits'salgorithm is evidently such an algorithm estimating Fourier components.Consequently, we have shown that the molecular vibronic spectra problemcorresponding to Fock-state boson sampling can be efficiently solved bya classical computer as accurately as running a boson sampler, whichindicates that there is no quantum advantage from this problem.

Gurvits' Algorithm with Repeated Rows and Columns

The disclosure below first presents the original Gurvits's algorithm andgeneralized algorithm to improve the bound when rows and columns arerepeated. The Gurvits's algorithm exploits the following equality:

$\begin{matrix}{{{Per}(X)} = {{{\mathbb{E}}_{x \in {\{{{- 1},1}\}}^{n}}\left\lbrack {\prod\limits_{i = 1}^{n}\left( {\sum\limits_{j = 1}^{n}{x_{i}X_{ij}x_{j}}} \right)} \right\rbrack}.}} & (37)\end{matrix}$

Thus, by sampling the random variable with uniform x∈{−1, 1}^(n), therandomized algorithm gives an estimate μ of permanent of an n×n matrix Xsuch that

|Per(X)−μ|<ϵ∥X∥ ^(n),  (38)

with high probability 1−δ with running time T=O(poly(n, 1/ϵ, log δ⁻¹)).

When a submatrix of a unitary matrix is considered, this expression issufficient because the submatrix's norm is bounded by 1. However, whenthe rows and columns in X from a unitary matrix are repeated, which isthe case when dealing with multiphoton inputs, the norm ∥X∥ is notgenerally bounded by 1. Therefore, the upper bound of the error canincrease as ∥X∥^(n), i.e., exponentially in n. To deal with suchmultiphoton cases with a reduced error, Gurvits's algorithm may need tobe generalized. Note that although the Gurvits's algorithm formultiphoton outputs may be generalized, i.e., Per(V_(n,m))=

n|{circumflex over (V)}|m

, where n∈{0, 1}^(n) and in is a nonnegative integer vector, that is notsufficient for the purpose of handling multiphoton input and output,i.e.,

n|{circumflex over (V)}|n

=Per (V_(n,n)).

Consider a case where the goal is to approximate a quantity

$\begin{matrix}{\frac{Per}{n!},} & (39)\end{matrix}$

where A∈

^(n×n) which is obtained by repeating ith row and column of a matrix B∈

^(k×k) for n_(i) times. Here, n∈

_(≥0) ^(k) with Σ_(i=1) ^(k) n_(i)=n. Further define a random variable

$\begin{matrix}{{{{{GenGly}_{x}(A)} \equiv {\frac{1}{\prod_{i = 1}^{k}n_{i}^{n_{i}}}{\prod\limits_{i = 1}^{k}{{\overset{\_}{y}}_{i}^{n_{i}}{\prod\limits_{i = 1}^{k}\left( {{y_{i}B_{i,1}} + \ldots + {y_{k}B_{i,k}}} \right)^{n_{i}}}}}}} = {\prod\limits_{i = 1}^{k}\left( \frac{{{\overset{\_}{y}}_{i}({By})}_{i}}{n_{i}} \right)^{n_{i}}}},} & (40)\end{matrix}$

where γ_(i)≡n_(i) ^(1/2)x_(i). Here, x∈

[n₁+1]× . . . ×

[n_(k)+1]≡

, where

[j] is the set of jth roots of unity. Then, its absolute value is givenas

$\begin{matrix}{{{❘{{GenGly}_{x}(A)}❘} = {{\left( {\prod\limits_{i = 1}^{k}n_{i}^{n_{i}}} \right)^{{- 1}/2}{❘{\prod\limits_{i = 1}^{k}\left( {{y_{i}b_{i,1}} + \ldots + {y_{k}b_{i,k}}} \right)^{n_{i}}}❘}} = {\left( {\prod\limits_{i = 1}^{k}n_{i}^{n_{i}}} \right)^{{- 1}/2}{❘{\prod\limits_{i = 1}^{k}c_{i}^{n_{i}}}❘}}}},} & (41)\end{matrix}$

where c_(i)≡γ₁B_(i,1)+ . . . +γ_(k)B_(i,k). Now, let us find its upperbound, which determines the error of the corresponding randomizedalgorithm by Chernoff bound. First, we have a constraint

$\begin{matrix}{{\sum\limits_{i = 1}^{k}{❘c_{i}❘}^{2}} = {{{\sum\limits_{i = 1}^{k}{❘({By})_{i}❘}^{2}} \leq {{B}^{2}{y}^{2}}} = {n{{B}^{2}.}}}} & (42)\end{matrix}$

Further considering:

$\begin{matrix}{{\log{\prod\limits_{i = 1}^{k}{❘c_{i}❘}^{n_{i}}}} = {\sum\limits_{i = 1}^{k}{n_{i}\log{{❘c_{i}❘}.}}}} & (43)\end{matrix}$

Using lagrange multiplier, with a constraint Σ_(i=1) ^(k)|c_(i)|²=c,

$\begin{matrix}{{Z \equiv {{\sum\limits_{i = 1}^{k}{n_{i}\log{❘c_{i}❘}}} - {\lambda\left( {{\sum\limits_{i = 1}^{k}{❘c_{i}❘}^{2}} - c} \right)}}},} & (44)\end{matrix}$

It may be obtained that:

$\begin{matrix}{{\frac{\partial Z}{\partial{❘c_{i}❘}} = {\frac{n_{i}}{❘c_{i}❘} = {{2\lambda{❘c_{i}❘}} = 0}}},{\forall{i.}}} & (45)\end{matrix}$

Using the constraint Eq. (42),

$\begin{matrix}{{{\sum\limits_{i = 1}^{k}{❘c_{i}❘}^{2}} = {{\sum\limits_{i = 1}^{k}\frac{n_{i}}{2\lambda}} = {\frac{n}{2\lambda} = c}}},} & (46)\end{matrix}$

the solution is obtained as

$\begin{matrix}{{❘c_{i}❘} = {\sqrt{\frac{n_{i}}{2\lambda}} = {\sqrt{\frac{{cn}_{i}}{n}}.}}} & (47)\end{matrix}$

Therefore, using the solution, we have an upper bound as

$\begin{matrix}{{{\prod\limits_{i = 1}^{k}{❘c_{i}❘}^{n_{i}}} \leq {\exp\left( {\sum\limits_{i = 1}^{k}{n_{i}\log\sqrt{\frac{{cn}_{i}}{n}}}} \right)}} = {{\prod\limits_{i = 1}^{k}{n_{i}^{n_{i}/2}\left( \frac{c}{n} \right)}^{n/2}} \leq {\prod\limits_{o = 1}^{k}{n_{i}^{n_{i}/2}{{B}^{n}.}}}}} & (48)\end{matrix}$

Hence, there is an upper bound of the absolute value of the randomvariable,

$\begin{matrix}{{❘{{GenGly}_{x}(A)}❘} = {{\left( {\prod\limits_{i = 1}^{k}n_{i}^{n_{i}}} \right)^{{- 1}/2}{❘{\prod\limits_{i = 1}^{k}c_{i}^{n_{i}}}❘}} \leq {{B}^{n}.}}} & (49)\end{matrix}$

Next, it may be proved that the average of the random variables gives anunbiased estimator of Per(A)/n!, i.e.,

$\begin{matrix}{{{{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {{GenGly}_{x}(A)} \right\rbrack} = {\frac{{Per}(A)}{n!}.{First}}},} & (50) \\{{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {\frac{1}{\prod\limits_{i = 1}^{k}n_{i}^{n_{i}}}{\prod\limits_{i = 1}^{k}{{\overset{\_}{y}}_{i}^{n_{i}}{\prod\limits_{i = 1}^{k}\left( {{y_{1}B_{i,1}} + \cdots + {y_{k}B_{i,k}}} \right)^{n_{i}}}}}} \right\rbrack} & (51) \\{= {\frac{1}{\prod\limits_{i = 1}^{k}n_{i}^{n_{i}}}{{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {\prod\limits_{i = 1}^{k}{{\overset{\_}{y}}_{i}^{n_{i}}\left( {y_{1},{B_{i,1} + \cdots + {y_{k}B_{i,k}}}} \right)}^{n_{i}}} \right\rbrack}}} & (52) \\{{= {\frac{1}{\prod\limits_{i = 1}^{k}n_{i}^{n_{i}}}{\sum\limits_{\sigma_{1},\ldots,{\sigma_{n} \in {\lbrack k\rbrack}}}{\left( {\prod\limits_{i = 1}^{k}B_{i,\sigma_{i}}} \right){{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {\prod\limits_{i = 1}^{k}{{\overset{\_}{y}}_{i}^{n_{i}}{\prod\limits_{i = 1}^{n}y_{\sigma_{i}}}}} \right\rbrack}}}}},} & (53)\end{matrix}$

where by the symmetry over the roots of unity,

$\begin{matrix}{{{{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {\prod\limits_{i = 1}^{k}{{\overset{\_}{y}}_{i}^{n_{i}}\prod\limits_{i = 1}^{n}}} \right\rbrack} = 0},} & (54)\end{matrix}$

unless the product insides is equal to Π_(i=1) ^(k)|γ_(i)|^(2n) ^(i) ,in which case it may be derived that

$\begin{matrix}{{{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {\prod\limits_{i = 1}^{k}{{\overset{\_}{y}}_{i}^{n_{i}}\prod\limits_{i = 1}^{n}}} \right\rbrack} = {\prod\limits_{i = 1}^{k}{n_{i}^{n_{i}}.}}} & (55)\end{matrix}$

Hence,

$\begin{matrix}{{{{\mathbb{E}}_{{\mathbb{x}} \in \chi}\left\lbrack {{GenGly}_{x}(A)} \right\rbrack} = {{\sum\limits_{{t_{1}\ldots},{{t_{k} \in {{\lbrack k\rbrack}:{❘{\{{{i:t_{i}} = j}\}}❘}}} = {n_{j}{\forall{j \in {\lbrack k\rbrack}}}}}}{\prod\limits_{i = 1}^{k}B_{i,t_{i}}}} = \frac{{Per}(A)}{n!}}},} & (56)\end{matrix}$

where n! appears to take into account the fact that each produce Π_(i=1)^(k)b_(i,t) _(k) appears n! times.

Therefore, GenGly_(x)(A) gives us an unbiased estimator of Per(A)/n!,and the random variables are bounded by ∥B∥^(n). Due to the Chernoffbound, the randomized algorithm with number of samples N=O(1/ϵ²)provides an estimate μ with high probability 1−δ such that

$\begin{matrix}{{❘{\frac{{Per}(A)}{n!} - \mu}❘} < {\epsilon{{B}^{n}.}}} & (57)\end{matrix}$

More Details of Exact Solution of Fourier Components of MolecularVibronic Spectra

A detailed derivation of the analytic solution of Fourier components ofthe vibronic spectra is described below.

The initial positive P-representation above may be represented by:

$\begin{matrix}{{P_{in}\left( {\alpha,\beta} \right)} = {\prod\limits_{i = 1}^{M}\left\lbrack {\frac{\sqrt{1 + \gamma_{i}}}{{\pi\gamma}_{i}}{\exp\left( {{{- \left( {\alpha_{i}^{2} + \beta_{i}^{2}} \right)}\left( {\gamma_{i}^{- 1} + {1/2}} \right)} + {\alpha_{i}\beta_{i}}} \right)}} \right\rbrack}} & (58) \\{{= {\mathcal{N}{\exp\left\lbrack {{{- \alpha} \cdot A \cdot \alpha} - {\beta \cdot B \cdot \beta} - {\alpha \cdot C \cdot \beta} - {\beta \cdot C^{T} \cdot \alpha}} \right\rbrack}}},} & (59)\end{matrix}$

where

$\begin{matrix}{{A \equiv {{diag}\left( {\gamma_{i}^{- 1} + {1/2}} \right)}_{i = 1}^{M}},{B \equiv {{diag}\left( {\gamma_{i}^{- 1} + {1/2}} \right)}_{i = 1}^{M}},{C \equiv {{- 1}/2}},{\mathcal{N} \equiv {\prod\limits_{i = 1}^{M}{\frac{\sqrt{1 + \gamma_{i}}}{{\pi\gamma}_{i}}.}}}} & (60)\end{matrix}$

Here, α and β are real M-dimensional vectors. In the disclosure above,the expression of the Fourier components of molecular vibronic spectrahas been found as:

{tilde over (G)}(k)=

e ^(−ikθ{circumflex over (n)}·w′)

,  (61)

where the expectation value is over the output state of Gaussian bosonsampling.

By using the relation e^(−γ{circumflex over (n)})=:e^({circumflex over (n)}(e) ^(−γ−1)) :, the operators can be written ina normal-ordered form and the Fourier transform may have the followingproperty:

$\begin{matrix}{{\overset{\sim}{G}(k)} = {\sum\limits_{\omega_{0} = 0}^{K}{\sum\limits_{m = 0}^{\infty}{{p(m)}{\delta\left( {\omega_{0} - {\omega \cdot m}} \right)}e^{{- {ik}}{\theta\omega}_{0}}}}}} & (62) \\{= {\sum\limits_{m = 0}^{\infty}{{p(m)}e^{{- {ik}}{{\theta\omega} \cdot m}}}}} & (63) \\\left. {{{= {\left\langle {\underset{i = 1}{\overset{M}{\otimes}}{\sum\limits_{m_{i} = 0}^{M}{❘m_{i}}}} \right\rangle\left\langle m_{i} \right.}}❘}e^{{- {ik}}{\theta\omega}_{i}m_{i}}} \right\rangle & (64) \\\left. {{{= {\left\langle {\underset{i = 1}{\overset{M}{\otimes}}{\sum\limits_{m_{i} = 0}^{M}{❘m_{i}}}} \right\rangle\left\langle m_{i} \right.}}❘}e^{{- {ik}}{\theta\omega}_{i}{\hat{n}}_{i}}} \right\rangle & (65) \\{= \left\langle {\underset{i = 1}{\overset{M}{\otimes}}e^{{- {ik}}{\theta\omega}_{i}{\hat{n}}_{i}}} \right\rangle} & (66) \\{{= \left\langle {:{\exp\left\lbrack {\sum\limits_{i = 1}^{M}{{\hat{n}}_{i}\left( {e^{{- {ikw}_{i}}\theta} - 1} \right)}} \right\rbrack}:} \right\rangle},} & (67)\end{matrix}$

where θ=2π/(K+1). For the last equality, we have used the relatione^(−γ{circumflex over (n)})=: e^({circumflex over (n)}(e) ^(−γ) ⁻¹⁾:.

The expectation value may be computed as:

$\begin{matrix}{{\overset{\sim}{G}(k)} = {\left\langle {:{\exp\left\lbrack {\sum\limits_{i = 1}^{M}{{\hat{n}}_{i}\left( {e^{i\phi_{i}} - 1} \right)}} \right\rbrack}:} \right\rangle = \left\langle {\exp\left\lbrack {\sum\limits_{i = 1}^{M}{\alpha_{i}^{\prime}{\beta_{i}^{\prime}\left( {e^{i\phi_{i}} - 1} \right)}}} \right\rbrack} \right\rangle_{P}}} & (68)\end{matrix}$ $\begin{matrix}{= \left\langle {\exp\left\lbrack {\sum\limits_{i = 1}^{M}{\left( {{\sum\limits_{j = 1}^{M}{U_{ij}\alpha_{j}}} + \alpha_{0i}} \right)\left( {{\sum\limits_{k = 1}^{M}{U_{ik}^{*}\beta_{k}}} + \beta_{0i}} \right)\left( {e^{i\phi_{i}} - 1} \right)}} \right\rbrack} \right\rangle_{P}} & (69)\end{matrix}$ $\begin{matrix}{= \left\langle {\exp\left\lbrack {{\alpha \cdot \mathcal{U} \cdot \beta} + {a \cdot \alpha} + {b \cdot \beta} + c_{0}} \right\rbrack} \right\rangle_{P}} & (70)\end{matrix}$ $\begin{matrix}{= {\mathcal{N}{\int{d\alpha d\beta{\exp\left\lbrack {{{- \alpha} \cdot A \cdot \alpha} - {\beta \cdot B \cdot \beta} - {\alpha \cdot \overset{\sim}{C} \cdot \beta} - {\beta \cdot {\overset{\sim}{C}}^{T} \cdot \alpha} + {a \cdot \alpha} + {b \cdot \beta} + c_{0}} \right\rbrack}}}}} & (71)\end{matrix}$ $\begin{matrix}{{= {\mathcal{N}{\int{{dq}{\exp\left\lbrack {{{{- q} \cdot Q \cdot q}/2} + {c \cdot q} + c_{0}} \right\rbrack}}}}},{{{where}\phi_{i}} \equiv {{- k}{\theta\omega}_{i}^{\prime}}},} & (72)\end{matrix}$ $\begin{matrix}{{\mathcal{U} \equiv {U^{T}\Phi U^{*}}},{\Phi \equiv {{diag}\left( {{e^{i\phi_{1}} - 1},\ldots,{e^{i\phi_{M}} - 1}} \right)}},{\overset{\sim}{C} \equiv {C - {\mathcal{U}/2}}},{Q = {2\begin{pmatrix}A & \overset{\sim}{C} \\{\overset{\sim}{C}}^{T} & B\end{pmatrix}}},} & (73)\end{matrix}$ $\begin{matrix}{{c \equiv \begin{pmatrix}a \\b\end{pmatrix}},{a = {U^{T}{\Phi\beta}_{0}}},{b = {U^{\dagger}{\Phi\alpha}_{0}}},{c_{0} = {\alpha_{0}{\Phi\beta}_{0}}},} & (74)\end{matrix}$

and

represents the average over the positive P-representation.

It is noted that the displacement occurs after the linear-opticalunitary, which matches the displacement δ/√{square root over (2)} in theDoktorov transformation, i.e., α₀=β₀=δ/√{square root over (2)}. It isfurther noted that Q's real part is positive-definite.

To see this, Q may be decomposed as:

Q = 2 ⁢ 2 ⊗ Γ - 1 + ( M - M - M M ) - ( 0 𝒰 𝒰 T 0 ) ( 75 ) = 2 ⁢ 2 ⊗ Γ -1 + ( M - U T ⁢ diag ⁡ ( e i ⁢ ϕ j ) j = 1 M ⁢ U * - U † ⁢ diag ⁡ ( e i ⁢ ϕ j )j = 1 M ⁢ U M ) , ( 76 )

where Γ≡diag(γ₁, . . . , γ_(M)).

The real part of Q may then be written as:

Re ⁡ ( Q ) = 2 ⁢ 2 ⊗ Γ - 1 + ( M L L T M ) , ( 77 )

where L≡Re(−U^(T) diag(e^(iϕj))_(j=1) ^(M)U*). Since ∥L∥_(op)≤1, Re(Q)is positive-definite, the integral always converges.

Using the Gaussian integral formula,

$\begin{matrix}{{{\int{{dq}{\exp\left\lbrack {{{{- q} \cdot Q \cdot q}/2} + {c \cdot q}} \right\rbrack}}} = {\sqrt{\frac{\left( {2\pi} \right)^{2M}}{❘Q❘}}{\exp\left( {\frac{1}{2}{c \cdot Q^{- 1} \cdot c}} \right)}}},} & (78)\end{matrix}$

the analytical form of the Fourier components may be obtained as:

G ~ ( k ) = 𝒩 ⁢ ( 2 ⁢ π ) 2 ⁢ M ❘ "\[LeftBracketingBar]" Q ❘"\[RightBracketingBar]" ⁢ exp ⁡ ( 1 2 ⁢ c · Q - 1 · c + c 0 ) , where ( 79) Q = ( 2 ⁢ Γ - 1 + M - U T ⁢ diag ⁡ ( e i ⁢ ϕ j ) j = 1 M ⁢ U * - U † ⁢ diag ⁡( e i ⁢ ϕ j ) j = 1 M ⁢ U 2 ⁢ Γ - 1 + M ) . ( 80 )

More Detail of the Positive P-Representation

The generalized (positive) P-representation introduced above isdescribed in further detail. Such generalized positive P representationmay be written as:

$\begin{matrix}{{\hat{\rho} = {\int_{{\mathbb{C}}^{2M}}{{P\left( {\alpha,\beta} \right)}{\hat{\Lambda}\left( {\alpha,\beta} \right)}d\alpha d\beta}}},{{\hat{\Lambda}\left( {\alpha,\beta} \right)} \equiv \frac{{\left. {❘\alpha} \right\rangle\left\langle \beta^{*} \right.}❘}{\left\langle {\beta^{*}❘\alpha} \right\rangle}},} & (81)\end{matrix}$

where |α

and |β*

are coherent states.

The generalized P-representation represents one of the quasi-probabilitydistributions of a bosonic state. An important property of thegeneralized P-representation is that the distribution P(α,β) can alwaysbe chosen non-negative (hence and name positive P-representation) andthe expectation value of a normal-ordered operator can be readilycomputed.

For simplicity, a simpler version may be used, where it is assumed thatthe Glauber-Sudarshan P-function exists, which is written as

{circumflex over (ρ)}=

d ² αP _(GS)(α)|α

α|,  (82)

for a single-mode state.

It is noted that although the existence of the Glauber-SudarshanP-function is assumed for simplicity, such assumption may be true forarbitrary density matrices. The assumption is that the positiveP-representation can be written as

$\begin{matrix}{{P\left( {\alpha,\beta} \right)} = {\frac{1}{4\pi^{2}}{\exp\left( {- \frac{{❘{\alpha - \beta^{*}}❘}^{2}}{4}} \right)}{\left\langle {\left( {\alpha + \beta^{*}} \right)/2{❘\hat{\rho}❘}\left( {\alpha + \beta^{*}} \right)/2} \right\rangle.}}} & (83)\end{matrix}$

If the assumption is true, then it guarantees that the representation ispositive and real.

By directly substituting the representation, it is obtained that:

$\begin{matrix}{{P\left( {\alpha,\beta} \right)} = {\frac{1}{4\pi^{2}}{\int{{P_{GS}(\gamma)}{\exp\left( {{{- \frac{1}{2}}{❘{\alpha - \gamma}❘}^{2}} - {\frac{1}{2}{❘{\beta^{*} - \gamma}❘}^{2}}} \right)}d^{2}{\gamma.}}}}} & (84)\end{matrix}$

In the meanwhile, for an analytical function ƒ:

$\begin{matrix}{{f(\gamma)} = {\frac{1}{2\pi}{\int_{\mathbb{C}}{{f(\alpha)}{\exp\left( {{- \frac{1}{2}}{❘{\alpha - \gamma}❘}^{2}} \right)}d^{2}{\alpha.}}}}} & (85)\end{matrix}$

It can be shown that, for example at γ=0,

$\begin{matrix}{{{\int_{\mathbb{C}}{{f(\alpha)}{\exp\left( {{- \frac{1}{2}}{❘\alpha ❘}^{2}} \right)}d^{2}\alpha}} = {{{2\pi{\exp\left( {\frac{1}{2}\frac{\partial^{2}}{{\partial\alpha}{\partial\beta^{*}}}} \right)}{f(\alpha)}}❘_{\alpha = 0}} = {2\pi{f(0)}}}},} & (86)\end{matrix}$

where the last equality holds because the function is analytic.

Using this property, the following may be obtained from Eq. (84)

$\begin{matrix}{\int_{{\mathbb{C}}^{2}}{{\hat{\Lambda}\left( {\alpha,\beta} \right)}{P\left( {\alpha,\beta} \right)}d^{2}\alpha d^{2}\beta}} & (87)\end{matrix}$ $\begin{matrix}{= {\frac{1}{4\pi^{2}}{\int_{{\mathbb{C}}^{3}}{d^{2}\alpha d^{2}\beta d^{2}\gamma{\hat{\Lambda}\left( {\alpha,\beta} \right)}{P_{GS}(\gamma)}{\exp\left( {{{- \frac{1}{2}}{❘{\alpha - \gamma}❘}^{2}} - {\frac{1}{2}{❘{\beta^{*} - \gamma}❘}^{2}}} \right)}}}}} & (88)\end{matrix}$ $\begin{matrix}{{{{\left. {= {\int_{\mathbb{C}}{d^{2}\alpha{P_{GS}(\alpha)}{❘\alpha}}}} \right\rangle\left\langle \alpha \right.}❘} = \hat{\rho}},} & (89)\end{matrix}$

which proves the assumption above.

It may be further shown that the expectation value of a normal-orderedoperator can be computed by averaging over the phase-space variable. Forexample, when computing for a single-mode: {circumflex over(α)}{circumflex over (α)}^(†){circumflex over (α)}{circumflex over(α)}^(†):=α^(†2){circumflex over (α)}²,

$\begin{matrix}{{{Tr}\left\lbrack {\hat{\rho}:\hat{a}{\hat{a}}^{\dagger}\hat{a}{\hat{a}}^{\dagger}:} \right\rbrack} = {{{Tr}\left\lbrack {\hat{\rho}{\hat{a}}^{\dagger 2}{\hat{a}}^{2}} \right\rbrack} = {\int_{{\mathbb{C}}^{2}}{{P\left( {\alpha,\beta} \right)}{{Tr}\left\lbrack {{\hat{\Lambda}\left( {\alpha,\beta} \right)}{\hat{a}}^{\dagger 2}{\hat{a}}^{2}} \right\rbrack}d\alpha d\beta}}}} & (90)\end{matrix}$ $\begin{matrix}{= {{\int_{{\mathbb{C}}^{2}}{{P\left( {\alpha,\beta} \right)}\frac{\left\langle {\beta^{*}{❘{{\hat{a}}^{\dagger 2}{\hat{a}}^{2}}❘}\alpha} \right\rangle}{\left\langle {\beta^{*}❘\alpha} \right\rangle}d\alpha d\beta}} = {\int_{{\mathbb{C}}^{2}}{{P\left( {\alpha,\beta} \right)}\beta^{2}\alpha^{2}d\beta d{\beta.}}}}} & (91)\end{matrix}$

Therefore, the expectation value may be directly integrated or estimatedby sampling (α,β) from P and averaging α²β² over the samples. The sameprocedure works for any normal-ordered operator ƒ({circumflex over (α)},{circumflex over (α)}^(†)) by replacing {circumflex over (α)}→α and{circumflex over (α)}^(†)→β. Especially when the operator is written asa function of {circumflex over (n)}, g({circumflex over (n)}), then thecorresponding phase-variable is given by g(αβ).

When applying a linear-optical circuit to a given input state, thequantum state transforms as

{circumflex over (ρ)}=

P(α,β){circumflex over (Λ)}(α,β)dαdβ→{circumflex over (ρ)}=

P(α,β){circumflex over (Λ)}(Uα,U*β)dαdβ.  (92)

which is written for the positive P-representation as

{circumflex over (ρ)}=

P(α,β){circumflex over (Λ)}(α,β)dαdβ→{circumflex over (ρ)}=

P(α,β){circumflex over (Λ)}(Uα,U*β)dαdβ.  (93)

Therefore, sampling from the positive P-representation of the outputstate of the linear-optical circuit is equivalent to sampling from thepositive P-representation of the input state and transform the phasevariable as (α,β)→(α′,β′)=(Uα,U*β).

Numerical Simulation

The various algorithms above may be numerically simulated to generatemolecular vibronic spectra.

The first example involves photo-electron spectra of formic acid (CH₂O₂,1¹A′→1²A′). To quantify the performance of the scheme based on Gaussianboson sampling and the proposed quantum-inspired classical algorithm,the ideal FCP distribution, generated by directly computing probabilitydistributions of the outcomes of Gaussian boson sampling, may becompared to the empirical distributions obtained by Gaussian bosonsampling simulation and the proposed classical algorithm.

It may be emphasized that the ideal Gaussian boson sampling assumes noerrors and losses. FIGS. 3 and 6 show that both schemes predict theground truth distribution well. In particular, FIG. 3 shows an exampledata and logic flow for the quantum-inspired classical algorithm of FIG.1 for generating molecular vibronic spectra, whereas FIG. 6 showscomparison of molecular vibronic spectra for CH₂O₂, 1¹A′→1²A′ asgenerated based on direct computation (ideal bars) and based on Gaussianboson sampling simulation (“GBS” points and curve).

In FIG. 7 , the total variation distances by varying the number ofsamples is shown. The total variance distance (TVD) is defined as

$\begin{matrix}{{TVD} \equiv {\frac{1}{2}{\sum\limits_{\omega_{vib}}{{❘{{{FCP}\left( \omega_{vib} \right)} - {{eFCP}\left( \omega_{vib} \right)}}❘}.}}}} & (94)\end{matrix}$

FIG. 7 particular shows comparison of total variance distance (TVD)between ideal direct computation (upper curve with circular points) andthe quantum-inspired classical algorithm (lower curve with triangularpoints) of FIG. 3 . FIG. 7 shows that while both schemes accuracyincreases as the number of samples increases, the proposed classicalalgorithm provides a better precision (about 5 times).

The second example is generated from a recent Gaussian boson samplingexperiments samples in order to verify that the quantum-inspiredclassical algorithm operates in the quantum supremacy regime.

Specifically, arbitrarily chosen virtual harmonic angular frequencies w′are introduced for the output modes to demonstrate that the algorithmcan even simulate a quantum supremacy regime for molecular vibronicspectra generation. For FIGS. 8 and 9 , nonzero weights are set only forthe first 6 modes to compare with its ideal distribution. More detailsfor this method are further disclosed below.

Since it corresponds to a Gaussian boson sampling on the small number ofmarginal modes, all ideal FC factors may be directly computed. FIGS. 8and 9 again show that both schemes predict the virtual spectracorrectly, while the spectra from the classical algorithm is moreaccurate (about 4 times) in terms of total variation distance.

Also, it is further observed that the total variance distance from theideal Gaussian boson sampling is smaller than the one from theexperimental Gaussian boson sampling, which is due to the experimentalnoise.

Finally, a random weight vector w′ has been generated for all 144 modesand compared the spectra from Gaussian boson sampling and classicalalgorithm, which is shown in FIG. 10 . While the spectra generated fromGaussian boson sampling fluctuates, which may be due to experimentalnoise, the classical algorithms spectra has the very similar shape toGaussian boson samplings spectra. Therefore, it confirms that thequantum-inspired classical algorithm gives a consistent result as theGaussian boson sampler in the quantum supremacy regime.

Approximation Algorithm of Hafnian

An example algorithm to approximate the hafnian above is provided.First, it is known that

$\begin{matrix}{{{x_{1}^{n_{1}}x_{2}^{n_{2}}\cdots x_{M}^{n_{M}}} = {\frac{1}{n!}{\sum\limits_{v_{1} = 0}^{n_{1}}{\cdots{\sum\limits_{v_{M} = 0}^{n_{M}}{\left( {- 1} \right)^{\sum\limits_{i = 1}^{M}v_{i}}\begin{pmatrix}n_{1} \\v_{1}\end{pmatrix}{\cdots\begin{pmatrix}n_{M} \\v_{M}\end{pmatrix}}\left( {\sum\limits_{i = 1}^{M}{h_{i}x_{i}}} \right)^{n}}}}}}},} & (95)\end{matrix}$

where h_(i)=n_(i)/2−v_(i).

Also a known proposition states that for z=(z₁, . . . , z_(M))^(T)˜

(0,Σ), where Σ is a real symmetric matrix and z can be a complexvariable vector,

$\begin{matrix}{{{{haf}\left( \Sigma_{n} \right)} = {{{\mathbb{E}}_{z}\left\lbrack {\prod\limits_{i = 1}^{M}z_{i}^{n_{i}}} \right\rbrack} = {\frac{1}{\left( {n/2} \right)!}{\sum\limits_{v_{1} = 0}^{n_{1}}{\cdots{\sum\limits_{v_{M} = 0}^{n_{M}}{\left( {- 1} \right)^{\sum\limits_{i = 1}^{M}v_{i}}\begin{pmatrix}n_{1} \\v_{1}\end{pmatrix}{\cdots\begin{pmatrix}n_{M} \\v_{M}\end{pmatrix}}\left( \frac{h^{T}\Sigma h}{2} \right)^{n/2}}}}}}}},} & (96)\end{matrix}$

where n=n₁+ . . . +n_(M) and h=(n₁/2−v₁, . . . , n_(M)/2−v_(M))^(T) andn is even.

It is noted that the equality between the hafnian and the sum in Eq.(96) holds even for general complex symmetric matrices. Especially whenn_(i)=1 for all i∈{1, . . . , n},

$\begin{matrix}{{{\mathbb{E}}_{z}\left\lbrack {\prod\limits_{i = 1}^{n}z_{i}} \right\rbrack} = {\frac{1}{\left( {n/2} \right)!}{\sum\limits_{v_{1} = 0}^{1}{\ldots{\sum\limits_{v_{M} = 0}^{1}{\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}\left( \frac{h^{T}{\sum h}}{2} \right)^{n/2}}}}}}} & (97)\end{matrix}$ $\begin{matrix}{= {\frac{1}{2^{n}}{\sum\limits_{v_{1} = 0}^{1}{\ldots{\sum\limits_{v_{n} = 0}^{1}{\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}\frac{2^{n}}{\left( {n/2} \right)!}\left( \frac{h^{T}{\sum h}}{2} \right)^{n/2}}}}}}} & (98)\end{matrix}$ $\begin{matrix}{{= {{\mathbb{E}}_{v}\left\lbrack {\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}\frac{2^{n/2}}{\left( {n/2} \right)!}\left( {h^{T}{\sum h}} \right)^{n/2}} \right\rbrack}},} & (99)\end{matrix}$

where h=(½−v₁, . . . , ½−v_(M))^(T). In this case, M=n.

Therefore,

$\begin{matrix}{\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}\frac{2^{n/2}}{\left( {n/2} \right)!}\left( {h^{T}{\sum h}} \right)^{n/2}} & (100)\end{matrix}$

is an unbiased estimator for the hafnian in the sense that

$\begin{matrix}{{{\mathbb{E}}_{v}\left\lbrack {\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}\frac{2^{n/2}}{\left( {n/2} \right)!}\left( {h^{T}{\sum h}} \right)^{n/2}} \right\rbrack} = {{{haf}(\sum)}.}} & (101)\end{matrix}$

Note that

$\begin{matrix}{{❘{\frac{2^{n/2}}{\left( {n/2} \right)!}\left( {h^{T}{\sum h}} \right)^{n/2}}❘} = {{❘{\frac{2^{n/2}}{\left( {n/2} \right)!}\left( {\frac{n}{4}{\overset{\hat{}}{h}}^{T}{\sum\overset{\hat{}}{h}}} \right)^{n/2}}❘} = {{❘{\frac{n^{n/2}}{\left( {n/2} \right){!2^{n/2}}}\left( {{\overset{\hat{}}{h}}^{T}{\sum\overset{\hat{}}{h}}} \right)^{n/2}}❘} \leq {\frac{\left( {n{\sum }} \right)^{n/2}}{\left( {n/2} \right){!2^{n/2}}}.}}}} & (102)\end{matrix}$

Therefore, in general, the estimate by sampling v and using theestimator has an error bound given by Chernoff bound as

$\begin{matrix}{{\Pr\left\lbrack {{❘{\overset{\sim}{p} - {{haf}(\sum)}}❘} > {\epsilon\frac{\left( {n{\sum }} \right)^{n/2}}{\left( {n/2} \right){!2^{n/2}}}}} \right\rbrack} \leq {2{{\exp\left( {{- \frac{1}{2}}\epsilon^{2}T} \right)}.}}} & (103)\end{matrix}$

Here, ∥Σ∥ is the spectral norm, i.e., the maximum singular value.

The Fourier components of molecular vibronic spectra, which is given by

$\begin{matrix}{{\overset{\sim}{G}(k)} = {\frac{{haf}\left( \sum_{n} \right)}{Z{n!}}.}} & (104)\end{matrix}$

Note that the operator norm of the matrix Σ_(α,β*) is 1 (see Eq. (??)).Since a matrix's submatrix's operator norm is smaller than or equal tothe matrix's norm, when the vector n is composed 0 or 1, Σ_(n)'soperator norm is smaller than or equal to 1.

In addition, when n has a component larger than 1, hafnian for a matrixobtained may be computed by repeating the row and column. Nevertheless,when repeating the row and column the operator norm does not increase asmuch as n!, i.e.,

$\begin{matrix}{\frac{\sum_{n}}{n!} \leq {{\sum }.}} & (105)\end{matrix}$

Therefore, the estimate by sampling x and using the estimator has anerror bound given by Chernoff bound as

$\begin{matrix}{{{\Pr\left\lbrack {{❘{\overset{\sim}{q} - \frac{{haf}(\sum)}{Z{n!}}}❘} \gtrsim {\epsilon\frac{e^{n/2}}{\sqrt{\pi n}Z}}} \right\rbrack} \leq {2{\exp\left( {{- \frac{1}{2}}\epsilon^{2}T} \right)}}},} & (106)\end{matrix}$

where Z=√{square root over (Π_(i=1) ^(n) cosh r_(i))} and we have usedStirling's formula, N!≈√{square root over (2πN)}(N/e)^(N) (moreprecisely, √{square root over (2πN)}(N/e)^(N)e^(1/12N+1)<N!<√{squareroot over (2πN)}(N/e)^(N)e^(1/12N)). Thus, if Z>e^(n/2), the estimationerror is smaller than e with exponentially small failure probability ifwe choose the number of samples as T=O(1/ϵ²).

For nonzero mean, one may again employ a similar formula for z˜N(μ, Σ):

$\begin{matrix}{{{lhaf}\left( \sum_{n} \right)} = {{\mathbb{E}}_{z}\left\lbrack {\prod\limits_{i = 1}^{M}z_{i}^{n_{i}}} \right\rbrack}} & (107)\end{matrix}$ $\begin{matrix}{{= {\sum\limits_{v_{1} = 0}^{n_{1}}{\ldots{\sum\limits_{v_{M} = 0}^{n_{M}}{\left( {- 1} \right)^{\sum_{i = 1}^{M}v_{i}}\begin{pmatrix}n_{1} \\v_{1}\end{pmatrix}\ldots\begin{pmatrix}n_{M} \\v_{M}\end{pmatrix}{\sum\limits_{r = 0}^{\lbrack{n/2}\rbrack}\frac{\left( \frac{h^{T}{\sum h}}{2} \right)^{r}\left( {h \cdot \mu} \right)^{n - {2r}}}{r{!{\left( {n - {2r}} \right)!}}}}}}}}},} & (108)\end{matrix}$

where n=n₁+ . . . +n_(M) and h=(n₁/2−v₁, . . . , n_(M)/2−v_(M))^(T) andΣ_(n) is obtained by replacing the diagonal elements of Σ by μ andrepeating ith row and column for n_(i) times. Again, let n_(i)=1 for alli's. In this case, M=n, and

$\begin{matrix}{{{lhaf}(\sum)} = {{{\mathbb{E}}_{z}\left\lbrack {\prod\limits_{i = 1}^{n}z_{i}} \right\rbrack} = {\sum\limits_{v_{1} = 0}^{1}{\ldots{\sum\limits_{v_{n} = 0}^{1}{\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}{\sum\limits_{r = 0}^{\lbrack{n/2}\rbrack}\frac{\left( \frac{h^{T}{\sum h}}{2} \right)^{r}\left( {h \cdot \mu} \right)^{n - {2r}}}{r{!{\left( {n - {2r}} \right)!}}}}}}}}}} & (109)\end{matrix}$ $\begin{matrix}{{= {\frac{1}{2^{n}\left\lbrack {n/2} \right\rbrack}{{{\mathbb{E}}_{v,r}\left\lbrack {\left( {- 1} \right)^{\sum_{i = 1}^{n}v_{i}}\frac{{2^{n}\left\lbrack {n/2} \right\rbrack}\left( \frac{h^{T}{\sum h}}{2} \right)^{r}\left( {h \cdot \mu} \right)^{n - {2r}}}{r{!{\left( {n - {2r}} \right)!}}}} \right\rbrack}.{Here}}}},} & (110)\end{matrix}$ $\begin{matrix}{{❘\frac{{2^{n}\left\lbrack {n/2} \right\rbrack}\left( \frac{h^{T}{\sum h}}{2} \right)^{r}\left( {h \cdot \mu} \right)^{n - {2r}}}{r{!{\left( {n - {2r}} \right)!}}}❘} = {\frac{\left\lbrack {n/2} \right\rbrack{n^{n/2}\left( {{\overset{\hat{}}{h}}^{T}{\sum\overset{\hat{}}{h}}} \right)}^{r}\left( {\overset{\hat{}}{h} \cdot \mu} \right)^{n - {2r}}}{2^{r}r{!{\left( {n - {2r}} \right)!}}} \leq \frac{\left\lbrack {n/2} \right\rbrack n^{n/2}{\Sigma }^{r}{\mu }^{n - {2r}}}{2^{r}r{!{\left( {n - {2r}} \right)!}}}}} & (111)\end{matrix}$

Therefore, the Chernoff bound provides

$\begin{matrix}{{\Pr\left\lbrack {{❘{p - {lh{{af}(\sum)}}}❘} > {\epsilon\frac{\left\lbrack {n/2} \right\rbrack n^{n/2}{\Sigma }^{r}{\mu }^{n - {2r}}}{2^{r}r{!{\left( {n - {2r}} \right)!}}}}} \right\rbrack} \leq {2{{\exp\left( {{- \frac{1}{2}}\epsilon^{2}T} \right)}.}}} & (112)\end{matrix}$

Especially when ∥Σ∥≤1, which is the case that we are interested in,

$\begin{matrix}{{\Pr\left\lbrack {{❘{p - {lh{{af}(\sum)}}}❘} > {\epsilon\frac{\left\lbrack {n/2} \right\rbrack n^{n/2}{\mu }^{n - {2r}}}{2^{r}r{!{\left( {n - {2r}} \right)!}}}}} \right\rbrack} \leq {2{{\exp\left( {{- \frac{1}{2}}\epsilon^{2}T} \right)}.}}} & (113)\end{matrix}$

The error bound clearly depends on the mean vector μ, which makes thebound more nontrivial than the hafnian case.

Other Techniques: Chernoff Bound

The accuracy of molecular vibronic spectra generation using Gaussianboson sampling may be derived. For example, let S={m_(i)}_(i=1) ^(M) bethe sample set, and for a fixed ω_(vib), the random variable may be theindicator I(m, ω_(vib)) i.e., S→I(S, ω_(vib))={I(m_(i),ω_(vib))}_(i=1)^(N), where

$\begin{matrix}{{I\left( {m,\omega_{vib}} \right)} = {{1{if}{w^{\prime} \cdot m}} \in \left\lbrack {{\omega_{vib} - \frac{{\Delta\omega}_{vib}}{2}},\ {\omega_{vib} + \frac{{\Delta\omega}_{vib}}{2}}} \right\rbrack}} & (114)\end{matrix}$ $\begin{matrix}{{I\left( {m,\omega_{vib}} \right)} = {0{{otherwise}.}}} & (115)\end{matrix}$

Thus, the Chernoff inequality for a given ω_(vib) may become

$\begin{matrix}{{{\Pr\left\lbrack {{❘{{{FCP}\left( \omega_{vib} \right)} - {\frac{1}{N}{\sum\limits_{m \in \mathcal{S}}{I\left( {m,\omega_{vib}} \right)}}}}❘} \geq \epsilon} \right\rbrack} \leq {2e^{{- 2}N\epsilon^{2}}}},} & (116)\end{matrix}$

where Σ_(m∈S)I(m, ω_(vib))/N is the estimate eFCP(ω_(vib)).

Then the Chernoff inequality may be expressed as:

Pr[|

( X )− X|≥ϵ]≤2e ^(−2Nϵ) ² ,  (117)

where X=(X₁+ . . . +X_(N))/N is the sample average and X_(i)∈[0, 1].

Methods for Numerical Results Above

The parameters for FIGS. 3 and 6-7 are provided in J. Huh, G. G.Guerreschi, B. Peropadre, J. R. McClean, and A. Aspuru-Guzik, Bosonsampling for molecular vibronic spectra, Nature Photonics 9, 615 (2015),which is herein incorporated by reference.

For FIGS. 8 and 9 above, the first 6 elements are selected as ω′ as[1000, 500, 300, 200, 400, 250] and zeros are used for the remainingelements. In addition, the raw data of the recent Gaussian bosonsampling is provided in https://quantum.ustc.edu.cn/web/node/951.

For the simulation, the frequencies in the unit of Δω_(vib) areapproximated depending on the number of bins and domain. As explainedabove, increasing the resolution requires an additional cost for thefast Fourier transform, which is only poly(K).

The generalized P-representation of a squeezed state is given by

$\begin{matrix}{{{P\left( {\alpha,\beta} \right)} = {\frac{\sqrt{1 + \gamma}}{\pi\gamma}e^{{{- {({\alpha^{2} + \beta^{2}})}}{({\gamma^{- 1} + {1/2}})}} + {\alpha\beta}}}},} & (118)\end{matrix}$

where α and β are real numbers, and γ=e^(2r) with r>0 being thesqueezing parameter.

The various classical algorithm above may be performed by any type ofcomputing system include a set of processors, memories, and othercomponents including but not limited to various user interfaces andnetwork units. At least one memory may computer instructions, whenexecuted by the set of processors, causes the computing system toperform the various algorithms described above.

Finally, it is to be understood that the following claims are intendedto cover all of the generic and specific features of the inventionherein described and all statements of the scope of the invention which,as a matter of language, might be said to fall therebetween.

What is claimed:
 1. A method for estimating a transition spectra in abosonic system having M bosonic modes, M being an integer equal to orlarger than 2, the method comprising: determining a unitary operationmodified from a transition operator associated with the bosonic system;processing each of N sampling bosonic states of the bosonic system in arepresentation space of the bosonic system to generate N sets of samplesusing at least the unitary operation, wherein each set of samplescorrespond to one of the N sampling bosonic states and are sampled fromat least one variable in the representation space, N being a positiveinteger; generating Fourier components of the transition spectra basedon the N sets of samples; and inverse-transforming the Fouriercomponents to estimate the transition spectra.
 2. The method of claim 1,wherein the unitary operation comprises a complex mixed position andmomentum operator.
 3. The method of claim 1, wherein each of the Nsampling bosonic states comprises a non-Gaussian state of the tM bosonicmodes.
 4. The method of claim 1, wherein each of the N sampling bosonicstates comprises a Gaussian state of the M bosonic modes.
 5. The methodof claim 1, wherein each of the N sampling bosonic states of the Mbosonic modes is expressed in a positive P-representation in therepresentation space of the bosonic system for the processing ofgenerating the corresponding set of samples in the representation space.6. The method of claim 5, wherein the M bosonic modes comprise Mmolecular vibronic modes.
 7. The method of claim 5, wherein the positiveP-representation of the each of the N sampling bosonic states of the Mbosonic modes and for the processing therein represents aquasi-probability distribution of the each of the N sampling bosonicstates of the M bosonic modes.
 8. The method of claim 7, wherein thephase space for the positive P-representation comprises 2M dimensions.9. The method of claim 7, wherein processing the each of the N samplingbosonic states of the M bosonic modes in the positive P-representationto generate the corresponding set of samples in the representation spacecomprises: selecting M pairs of complex numbers for two representationspace variables of the each of the N sampling bosonic states in therepresentation space; transforming the M pairs of complex numbers byapplying a unitary matrix corresponding to the unitary operation togenerate a transformed M pairs of complex numbers; and generating thecorresponding set of samples from the transformed M pairs of complexnumbers.
 10. The method of claim 9, wherein the transition operatorassociated with the bosonic system is decomposed into the unitaryoperation, a dressing operation, and a displacement operation.
 11. Themethod of claim 9, wherein at least one of the two phase space variablesis associated with mean quantum numbers of M-mode coherent states of thebosonic system.
 12. The method of claim 1, wherein: the transitionoperator associated with the bosonic system is decomposed into theunitary operation, a dressing operation, a displacement operation, and atwo-mode squeezing operation; M comprises a even integer; the bosonicsystem comprises a molecular vibronic system; and the M bosonic modescomprises M/2 vibronic modes of the molecular vibronic system and M/2auxiliary modes.
 13. A computing system for estimating a transitionspectra in a bosonic system having M bosonic modes, M being an integerequal to or larger than 2, the system comprising one or more processorsand a memory for storing computer instructions, the one or moreprocessors, when executing the computer instructions, are configured tocause the computing system to: determine a unitary operation modifiedfrom a transition operator associated with the bosonic system; processeach of N sampling bosonic states of the bosonic system in arepresentation space of the bosonic system to generate N sets of samplesusing at least the unitary operation, wherein each set of samplescorrespond to one of the N sampling bosonic states and are sampled fromat least one variable in the representation space, N being a positiveinteger; analytically generate Fourier components the transition spectrabased on the N sets of samples; and inverse-transform the Fouriercomponents to estimate the transition spectra.
 14. The computing systemof claim 13, wherein the unitary operation comprises a complex mixedposition and momentum operator.
 15. The computing system of claim 13,wherein each of the N sampling bosonic states comprises a non-Gaussianstate of the tM bosonic modes.
 16. The computing system of claim 13,wherein each of the N sampling bosonic states comprises a Gaussian stateof the M bosonic modes.
 17. The computing system of claim 16, whereineach of the N sampling bosonic states of the M bosonic modes isexpressed in a positive P-representation in the representation space ofthe bosonic system for the processing of generate the corresponding setof samples in the representation space.
 18. The computing system ofclaim 17, wherein: the M bosonic modes comprise M molecular vibronicmodes, and wherein the positive P-representation of the each of the Nsampling bosonic states of the M bosonic modes and for the processingtherein represents a quasi-probability distribution of the each of the Nsampling bosonic states of the M bosonic modes; and the phase space forthe positive P-representation comprises 2M dimensions.
 19. The computingsystem of claim 18, wherein to process the each of the N samplingbosonic states of the M bosonic modes in the positive P-representationto generate the corresponding set of samples in the representation spacecomprises: select M pairs of complex numbers for two representationspace variables of the Gaussian state in the representation space;transform the M pairs of complex numbers by applying s unitary matrixcorresponding to the unitary operation to generate a transformed M pairsof complex numbers; and generate the corresponding set of samples fromthe transformed M pairs of complex numbers.
 20. The computing system ofclaim 19, wherein: the transition operator associated with the bosonicsystem is decomposed into the unitary operation, a dressing operation,and a displacement operation; and at least one of the two phase spacevariables is associated with mean quantum numbers of M-mode coherentstates of the bosonic system.